{ "cells": [ { "cell_type": "markdown", "id": "7aaaac4f", "metadata": {}, "source": [ "# Copula estimation with fastKDE\n", "\n", "```{warning}\n", "This capability is experimental and relatively untested.\n", "```\n", "\n", "This notebook demonstrates using `fastKDE` to estimate the copula density of a given dataset." ] }, { "cell_type": "code", "execution_count": 1, "id": "7c91dae7", "metadata": {}, "outputs": [], "source": [ "\"\"\" Import libaries \"\"\"\n", "\n", "import fastkde.fastKDE as fastKDE\n", "import fastkde.plot\n", "import numpy as np\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "8c2b2227", "metadata": {}, "outputs": [], "source": [ "\"\"\" Give a simple example of calculating a copula. \"\"\"\n", "\n", "# Generate some data\n", "np.random.seed(0)\n", "n = int(1e6)\n", "\n", "# sample from a bivariate normal distribution with zero correlation;\n", "# the copula PDF of this distribution is uniform\n", "xy = np.random.multivariate_normal([0, 0], [[1, 0.0], [0.0, 1]], n)\n", "x = xy[:, 0]\n", "y = xy[:, 1]\n", "\n", "# initialize the fastKDE object\n", "kde = fastKDE.fastKDE([x,y])\n", "\n", "# calculate the copula\n", "copula = kde.getCopula()" ] }, { "cell_type": "code", "execution_count": 3, "id": "2f83c654", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAH5CAYAAAB9Dt2yAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAi6dJREFUeJzt3Xl4U8X6B/DvyV5KW5ZCF2lLQWUXsSAFRMWlLIriCqgIV+XKxQ24/LzickX0igsqoiyCLOIGXhFXlEXZFJBF8IogixRaoLW2QguUpm0yvz/O0uYkaRvoSZry/TzPPOk5mTOZJJBO37wzIwkhBIiIiIiIwowp1B0gIiIiIjoTHMgSERERUVjiQJaIiIiIwhIHskREREQUljiQJSIiIqKwxIEsEREREYUlDmSJiIiIKCxZQt2B2uZ2u3H06FFERUVBkqRQd4eIiIhqQAiBEydOIDExESZT8ONsJSUlKC0tNax9m80Gh8NhWPvnqno3kD169CiSkpJC3Q0iIiI6A9nZ2WjRokVQH7OkpASpKQ2Rm+cy7DHi4+ORmZnJwWwtq3cD2aioKADyf4To6OgQ94aIiIhqoqioCElJSdrv8WAqLS1Fbp4Lh7a1RHRU7UeDi064kZJ2EKWlpRzI1rJ6N5BV0wmio6M5kCUiIgozoUwLbBgloWFU7T++G0x1NAonexERERFRWKp3EVkiIiKiM+ESbriEMe2SMRiRJSIiIqKwxIgsEREREQA3BNyo/ZCsEW2SjANZIiIiIgBuuGFEEoAxrRLA1AIiIiIiClOMyBIREREBcAkBl6j9NAAj2iQZI7JEREREdci6deswcOBAJCYmQpIkfPrpp1XWz8nJwR133IE2bdrAZDJhzJgxXnXmzJmD3r17o3HjxmjcuDGuueYabN682aPOxIkTIUmSR4mPj6/FZ1b7OJAlIiIiQsVkLyNKIE6dOoXOnTvjzTffrFF9p9OJZs2a4YknnkDnzp191lmzZg2GDh2K1atXY+PGjUhOTkZGRgaOHDniUa9Dhw7IycnRyi+//BJQ34ONqQVEREREdUj//v3Rv3//Gtdv2bIlXn/9dQDAvHnzfNZ5//33PY7nzJmDjz/+GN9++y3uvvtu7bzFYqnzUdjKGJElIiIighyRdRlQ1IhsUVGRR3E6nSF7rsXFxSgrK0OTJk08zu/btw+JiYlITU3FkCFDcODAgRD1sGY4kCUiIiIKgqSkJMTExGhl8uTJIevLY489hvPOOw/XXHONdq579+5YuHAhli9fjjlz5iA3Nxc9e/ZEQUFByPpZHaYWEBEREcH4DRGys7MRHR2tnbfb7bX+WDXx0ksv4cMPP8SaNWvgcDi085XTGTp16oQePXqgdevWeOeddzBu3LhQdLVaHMgSERERBUF0dLTHQDYUpkyZgueffx6rVq3CRRddVGXdyMhIdOrUCfv27QtS7wLH1AIiIiIiVKwja0SpC15++WU8++yz+Oabb9C1a9dq6zudTuzevRsJCQlB6N2ZYUSWiIiICIBbKUa0G4iTJ09i//792nFmZiZ27NiBJk2aIDk5GRMmTMCRI0ewcOFCrc6OHTu0a//880/s2LEDNpsN7du3ByCnEzz11FP44IMP0LJlS+Tm5gIAGjZsiIYNGwIAxo8fj4EDByI5ORl5eXl47rnnUFRUhOHDh5/5kzcYB7JEREREdcjWrVvRp08f7VjNTx0+fDgWLFiAnJwcZGVleVzTpUsX7edt27bhgw8+QEpKCg4ePAgAmDFjBkpLS3Hrrbd6XPf0009j4sSJAIDDhw9j6NChyM/PR7NmzZCeno5NmzYhJSXFgGdZOyQh6ki8u5YUFRUhJiYGhYWFIc9DISIiopoJ5e9v9bF/3d0cUVG1n3V54oQbHdrlcWxiAObIEhEREVFYYmoBEREREQCXkIsR7ZIxGJElIiIiorDEiCwRERER6s6qBVRzjMgSERERUVhiRJaIiIgIgBsSXJAMaZeMwYEsEREREQC3kIsR7ZIxmFpARERERGGJEVkiIiIiAC6DUguMaJNkjMgSERERUVhiRJaIiIgIjMiGI0ZkiYiIiCgsMSJLREREBMAtJLiFActvGdAmyRiRJSIiIqKwxIgsEREREZgjG444kCUiIiIC4IIJLgO+rHbVeoukYmoBEREREYUlRmSJiIiIAAiDJnsJTvYyDCOyRERERBSWGJElIiIiAid7hSNGZImIiIgoLDEiS0RERATAJUxwCQNWLRC13iQpGJElIiIiorDEiCwRERERADckuA2I8bnBkKxROJAlIiIiAid7hSOmFhARERFRWGJEloiIiAhGTvZiaoFRGJElIiIiorBk6EC2ZcuWkCTJqzzwwAM+669Zs8Zn/d9++83IbhIREREpk72MKWQMQ1MLtmzZApfLpR3v3LkT1157LW677bYqr9uzZw+io6O142bNmhnWRyIiIiIKT4YOZPUD0BdeeAGtW7fGFVdcUeV1zZs3R6NGjQzsGREREZEnN0xwcfmtsBK0HNnS0lK89957uOeeeyBJVYfYu3TpgoSEBFx99dVYvXp1lXWdTieKioo8ChERERHVf0EbyH766ac4fvw4RowY4bdOQkICZs+ejSVLluCTTz5BmzZtcPXVV2PdunV+r5k8eTJiYmK0kpSUZEDviYiIqL5TVy0wopAxJCGCsyZE3759YbPZ8MUXXwR03cCBAyFJEj7//HOf9zudTjidTu24qKgISUlJKCws9MizJSIiorqrqKgIMTExIfn9rT72Bzs6okGUudbbLz7hwh0X7+TYxABBWUf20KFDWLVqFT755JOAr01PT8d7773n93673Q673X423SMiIiKiMBSUgez8+fPRvHlzXHfddQFfu337diQkJBjQKyIiIqIKLiHBJQzYotaANklm+EDW7XZj/vz5GD58OCwWz4ebMGECjhw5goULFwIApk6dipYtW6JDhw7a5LAlS5ZgyZIlRneTiIiIiMKM4QPZVatWISsrC/fcc4/XfTk5OcjKytKOS0tLMX78eBw5cgQRERHo0KEDvvrqKwwYMMDobhIREdE5zmXQ8lsuLr9lmKBN9gqWUCaLExER0ZmpC5O9FmzvbNhkrxFdfubYxABByZElIiIiquvcwgS3AUtluetXzLBO4cJmRERERBSWGJElIiIiAnNkwxEHskREREQA3DBmqSx3rbdIKqYWEBEREVFYYkSWiIiICIAbJrgNiPEZ0SbJ+MoSERERUVhiRJaIiIgIgEuY4DJg+S0j2iQZX1kiIiIiCkuMyBIREREBcEOCG0asWlD7bZKMEVkiIiIiCkuMyBIRERGBObLhiK8sERERESp29jKiBGLdunUYOHAgEhMTIUkSPv300yrr5+Tk4I477kCbNm1gMpkwZswYn/WWLFmC9u3bw263o3379li6dKlXnRkzZiA1NRUOhwNpaWlYv359QH0PNg5kiYiIiOqQU6dOoXPnznjzzTdrVN/pdKJZs2Z44okn0LlzZ591Nm7ciMGDB2PYsGH4+eefMWzYMNx+++348ccftTqLFy/GmDFj8MQTT2D79u3o3bs3+vfvj6ysrFp5XkaQhBD1agPgoqIixMTEoLCwENHR0aHuDhEREdVAKH9/q4/90pbeiGhY+1mXp0+W49Fu68/ouUmShKVLl2LQoEE1qn/llVfi4osvxtSpUz3ODx48GEVFRfj666+1c/369UPjxo3x4YcfAgC6d++OSy65BDNnztTqtGvXDoMGDcLkyZMD6newMCJLREREFARFRUUexel0Bu2xN27ciIyMDI9zffv2xYYNGwAApaWl2LZtm1edjIwMrU5dxIEsEREREeStZI3Ij1W3qE1KSkJMTIxWghnlzM3NRVxcnMe5uLg45ObmAgDy8/PhcrmqrFMXcdUCIiIioiDIzs72SC2w2+1BfXxJ8lzPVgjhda4mdeoSDmSJiIiIALiFCW4DlspS24yOjg7Z/J34+HivyGpeXp4WgY2NjYXZbK6yTl3E1AIiIiKieq5Hjx5YuXKlx7kVK1agZ8+eAACbzYa0tDSvOitXrtTq1EWMyBIREREBcEGCy4DtZANt8+TJk9i/f792nJmZiR07dqBJkyZITk7GhAkTcOTIESxcuFCrs2PHDu3aP//8Ezt27IDNZkP79u0BAI888gguv/xyvPjii7jxxhvx2WefYdWqVfj++++1NsaNG4dhw4aha9eu6NGjB2bPno2srCyMGjXqLJ69sTiQJSIiIoLxqQU1tXXrVvTp00c7HjduHABg+PDhWLBgAXJycrzWdu3SpYv287Zt2/DBBx8gJSUFBw8eBAD07NkTixYtwpNPPomnnnoKrVu3xuLFi9G9e3ftusGDB6OgoACTJk1CTk4OOnbsiGXLliElJSXQpxw0XEeWiIiIQq4urCP7zI/XwGHAOrIlJ8vxdPdVHJsYgBFZIiIiIgAuBJ4GUNN2yRic7EVEREREYYkRWSIiIiLUnRxZqjm+skREREQUlhiRJSIiIgLgEia4DIieGtEmyfjKEhEREVFYYkSWiIiICICABLcBqxYIA9okGSOyRERERBSWGJElIiIiAnNkwxEHskREREQA3EKCW9R+GoARbZKMfyIQERERUVhiRJaIiIgIgAsmuAyI8RnRJsn4yhIRERFRWGJEloiIiAjMkQ1HjMgSERERUVhiRJaIiIgIgBsmuA2I8RnRJsn4yhIRERFRWGJEloiIiAiAS0hwGZDPakSbJONAloiIiAic7BWOmFpARERERGGJEVkiIiIiAEKY4Ba1H+MTBrRJMr6yRERERBSWGJElIiIiAuCCBBcMmOxlQJskY0SWiIiIiMISI7JEREREANzCmBUG3KLWmySFoRHZiRMnQpIkjxIfH1/lNWvXrkVaWhocDgdatWqFWbNmGdlFIiIiIgpThkdkO3TogFWrVmnHZrPZb93MzEwMGDAAI0eOxHvvvYcffvgBo0ePRrNmzXDLLbcY3VUiIiI6h7kNWrXAiDZJZvhA1mKxVBuFVc2aNQvJycmYOnUqAKBdu3bYunUrpkyZ4ncg63Q64XQ6teOioqKz7jMRERGde9yQ4DZgYpYRbZLM8D8R9u3bh8TERKSmpmLIkCE4cOCA37obN25ERkaGx7m+ffti69atKCsr83nN5MmTERMTo5WkpKRa7T8RERER1U2GDmS7d++OhQsXYvny5ZgzZw5yc3PRs2dPFBQU+Kyfm5uLuLg4j3NxcXEoLy9Hfn6+z2smTJiAwsJCrWRnZ9f68yAiIqL6zyUkwwoZw9DUgv79+2s/d+rUCT169EDr1q3xzjvvYNy4cT6vkSTPN1sI4fO8ym63w26311KPiYiIiChcBHX5rcjISHTq1An79u3zeX98fDxyc3M9zuXl5cFisaBp06bB6CIRERGdozjZK/wE9ZV1Op3YvXs3EhISfN7fo0cPrFy50uPcihUr0LVrV1it1mB0kYiIiIjChKED2fHjx2Pt2rXIzMzEjz/+iFtvvRVFRUUYPnw4ADm/9e6779bqjxo1CocOHcK4ceOwe/duzJs3D3PnzsX48eON7CYRERGRvGqBMKBw1QLDGJpacPjwYQwdOhT5+flo1qwZ0tPTsWnTJqSkpAAAcnJykJWVpdVPTU3FsmXLMHbsWEyfPh2JiYmYNm0a15AlIiIiIi+GDmQXLVpU5f0LFizwOnfFFVfgp59+MqhHRERERL4Jg9aRFYzIGiaok72IiIiI6io1FcCIdskYnEZHRERERGGJEVkiIiIicPmtcMRXloiIiIjCEiOyRERERGCObDhiRJaIiIiIwhIjskRERERQNkQwYKksbohgHEZkiYiIiCgsMSJLREREBObIhiMOZImIiIjAgWw4YmoBEREREYUlRmSJiIiIwIhsOGJEloiIiIjCEiOyRERERGBENhwxIktERERUh6xbtw4DBw5EYmIiJEnCp59+Wu01a9euRVpaGhwOB1q1aoVZs2Z53H/llVdCkiSvct1112l1Jk6c6HV/fHx8bT+9WsWILBEREREAAWM2LxAB1j916hQ6d+6Mv/3tb7jllluqrZ+ZmYkBAwZg5MiReO+99/DDDz9g9OjRaNasmXb9J598gtLSUu2agoICdO7cGbfddptHWx06dMCqVau0Y7PZHGDvg4sDWSIiIqI6pH///ujfv3+N68+aNQvJycmYOnUqAKBdu3bYunUrpkyZog1kmzRp4nHNokWL0KBBA6+BrMViqfNR2MqYWkBERESEihxZIwoAFBUVeRSn01kr/d64cSMyMjI8zvXt2xdbt25FWVmZz2vmzp2LIUOGIDIy0uP8vn37kJiYiNTUVAwZMgQHDhyolT4ahQNZIiIiIhg/kE1KSkJMTIxWJk+eXCv9zs3NRVxcnMe5uLg4lJeXIz8/36v+5s2bsXPnTtx3330e57t3746FCxdi+fLlmDNnDnJzc9GzZ08UFBTUSj+NwNQCIiIioiDIzs5GdHS0dmy322utbUnyzO0VQvg8D8jR2I4dO+LSSy/1OF85naFTp07o0aMHWrdujXfeeQfjxo2rtb7WJg5kiYiIiGD88lvR0dEeA9naEh8fj9zcXI9zeXl5sFgsaNq0qcf54uJiLFq0CJMmTaq23cjISHTq1An79u2r1f7WJqYWEBEREYWxHj16YOXKlR7nVqxYga5du8JqtXqc/+ijj+B0OnHXXXdV267T6cTu3buRkJBQq/2tTRzIEhEREcH4HNmaOnnyJHbs2IEdO3YAkJfX2rFjB7KysgAAEyZMwN13363VHzVqFA4dOoRx48Zh9+7dmDdvHubOnYvx48d7tT137lwMGjTIK1ILAOPHj8fatWuRmZmJH3/8EbfeeiuKioowfPjwgPofTEwtICIiIqpDtm7dij59+mjHan7q8OHDsWDBAuTk5GiDWgBITU3FsmXLMHbsWEyfPh2JiYmYNm2a1xq0e/fuxffff48VK1b4fNzDhw9j6NChyM/PR7NmzZCeno5NmzYhJSXFgGdZOyShZgPXE0VFRYiJiUFhYaEheShERERU+0L5+1t97F6fPQhLZO1NwFKVn3Lihxvf5NjEAEwtICIiIqKwxNQCIiIiIsjb0xqxRa0RbZKMA1kiIiIiGL/8FtU+phYQERERUVhiRJaIiIgIgBAShAHRUyPaJBkjskREREQUlhiRJaK6qbQU+OsvuRQVyeXUKaCkBCgrA1wuwGQCzGbAbgciIoDoaLk0agQ0bQo0bAj42GeciMgX5siGHw5kiSg0XC4gMxP49Vdg715g/37g0CEgOxvIyQGOHTv7x3A4gPh4oEULIDkZaN0auPBCoH17oF07efBLRERhiwNZIjKeEPIAdfNmYONG+Xb7djnCaqSSEuDgQbnomc3yYLZbN6BnT+DSS4EOHeTzRHROYo5s+OFAloiMkZMDLF8OrFoFrF0LHD5co8tKrBbkxsQgP7ohjkVGorBhBE5EOFAU6UCxzYZSuxllFjPcJgkSALPLBXt5OSJKy9DwtBNRxSWIKi5BkxOn0LToFJofP4HGp4q9H8jlAnbulMv8+fK5qCigVy/gqquAfv2Ajh2ZmkBEVIdxIEtEtUMIYNcu4IsvgE8/BX78scrqR2IaY3f8efgtIQEHYpvjYHwsDjdpgvzGkYAkwW11yxUt8i7aQrmVzAKQlJ21lTGmJHnutK1FP5RbR3EZWuQfQ0ruX2idm4cLj/yB9oeP4oKjf8DqdldceOIE8M03cnn0USAlBbjxRuD664ErrgBstjN/fYiozhMG5cgyImscDmSJ6Oz8/jvw7rvAokXAnj0+qxRbbPileQp+aZ6C/yWm4H/xKSho1BAA4FK2NXcrY0SLkm3gtspf8asDWKF84y9MAkJdb8Wk3Oevb275l4fTbcHvMRE4EHUeVl8ASOXy3Y7TZeh0+DAuPpSFztmHkHYoE81PFlVcf+gQMG2aXBo1kge1d90lR2xNXPSFqL4RkP8mN6JdMgYHskQUuFOngI8/Bt55B1i92meV36PjsT6+AzYmt8XO2GQ4I+WPm3JH3YlMOK02bE1thZ9atAIAmEoFUgvy0HvvHlyxfxcuzfodVrdLrnz8uPx833kHSEoChg2TB7Xt2oXuCRARneM4kCWimtu/X45OLlggfw1fiRvAzoYpWB/XERubtsXhJs0BAOURJkhOwKJETyXlm3yTS1JuletL5VuXEpkVyqeT2yLX0yKyZqkiIqukFPj71k7LOFAis+pjS8pj6o9NLgmHGsUhu0scPuhyOaJPnsZlB/agz4Ff0PvgLjQsc8oVs7OB55+Xy2WXAf/4B3DrrUw9IApzbkiQYMDyWwa0STIOZImoakLIE7beeAPuL77w2kXlsLUxVsRchFXNu6DAFg13hBUAYHbKo8OKQaZ8pZbeqg0i5Qpuq1JL+drfbfG8rRjIotJAVvJ8DP3vCvWxtAGt57F+IKvdr9yedERgeduLsbztxbA7y3D5wV24cddm9MjeA4tQKn3/vVzGjwf+/nd5UBsXByIiMh4HskTkm9sNfPYZ8NxzwE8/AajYCrAEZqyxn4/lMZ3xW0xLeXKWzR6yrgZDqcWKVed3xnctO6NJ8Qn02/sTbvptE1of/0OukJMDPPMMMHmynHbw2GPA+eeHttNEFBAuvxV+OJAlIk8uF/Dhh8ALL8ibFVTypxSJL+zt8HX0RThpcgB2mxYVlZTZ/1K5eiufN5epSwuorSiRWd3X+y71WInQSsqnk0mNzJoqLnerS73qIrHVphjobvWRWf39Hm0qj+kyAX/aovDupVfg3W6XI+3I77jz5/W44sCvcpS2tBSYO1de0uvWW4GnnpKX8SIiolrHgSwRyYSQl856/HGvAexeqQkWWzpgY4Pz4ZZMkEz1O/paY5KEbS3Ox7YW56P5ieO44+f1uHXnRkSVlsgR7Y8+ksvttwOTJgFt2oS6x0RUBbeQIHGL2rDCgSwRAWvWABMmAJs2eZz+1RyHRY7O2OpoKUdeLWY5+GlRw6WVdsFSI7JKaNWkRGKFSZnUpW0soIY/PbNtJWVClsutTApTIrNCyZlVI7TCVLHylZorqwVQq/tdoVsDR/J33kdE1uuxdH042qARpmQMxMw+1+CObT/gzi3r0eyUMiHuo4/kVR6GD5dTNRITq+koERHVBBdCJDqX7d8P3HQT0KePxyB2t6kZHo3oi39GXY8t1iTubhWAk44IzO51Da59+Ak8n3Ej8iPl9XLhdsvpBhdeKA9mT58ObUeJyIsQxhUyBiOyROeiU6eA//wHeOUVOadTcdDSBAtjemCjoxUg+ViERh3QmpW/gU0+BrhqwNUlf3KbypVPcJN8h1DCqeYyf5/sSo6s8snv1lYSqFiGS+gjsv5WLahOdb9cpIomtcdQHlPN0xXKagvqqgvqKgulNivmXXM5PrjmUgxf9z3+/u0axJw+Lb/2Tz0FzJkDTJki59HyDwUiojPCiCzRueaLL4D27eXZ9cog9i9TA7zS5Fo8EHcHNka05sCqFhXb7Zh57dW4YuK/sOCKXihX8yKysuTc2YEDgYMHQ9pHIpKpqxYYUcgYhg5kJ0+ejG7duiEqKgrNmzfHoEGDsMfPFpaqNWvWQJIkr/Lbb78Z2VWi+i8nBxg0CLjhBnkQBaBMMmFR/GW4t/2DWJmUDlfjaIhGDYGGDeQS4ZCL3SYXm1UuFotcTCY5KlspMisJ4VHg9iwml1ygFPXYVK6UMqWUQ1fk+6VyaMW7jmeRXL6LVsflWSS3d5H3rPTxekpyESbP4rbKxeUQcDkE3A433A43jsfbMfGe6zFwygNY17nSslxffQV06gS89Ra/fyQKMQ5kw4+hqQVr167FAw88gG7duqG8vBxPPPEEMjIysGvXLkRGRlZ57Z49exAdHa0dN2vWzMiuEtVfQgDvvQc88ghw7Jh2emvT8zHtghtwuEEzmMqUCVql8q2kbGagTtxSl9TSBlpuNV1AqojeapsT+P7AVidQCS1VQG1Lt+uW7n5Ja195SAitG9ovB31qgX7FrxrSpxEIU6WluapZ4qviGuFZ3yIfW2xyQ0fOb4QHJw3BDT/9D2Pe/BaxBaeAkyeBUaOApUuBefM4GYyIqIYMHch+8803Hsfz589H8+bNsW3bNlx++eVVXtu8eXM0atTIwN4RnQP++kseIP33v9qpAntDvNHuBqyOvwjmUkYAQ0KSsK73hdh6SQrGvv0t+n6xSz6/fLkcnZ07V46eE1FQcfmt8BPUHNnCwkIAQJMmTaqt26VLFyQkJODqq6/G6tWr/dZzOp0oKiryKEQEYN065Ca39hjELk/tgiGDHsWKCy9BWbQF5Q3NKG9ohivCAleEBW67XIRDKVazZ7GYlGKWi8kEYZYgzEpkVpLkTxWTHJn1F52tkjbNF4CQI7mS8PFVv5CjtpJbaHW0+9y6InRFf7+uXuXH1F+j749WtPOSPDFNn5KgVDCZ5GI1uVAWZcF7T/TAS6/1RVEzh1zvr7/klST++U+grCzw14+I6BwStIGsEALjxo3DZZddho5V7HKTkJCA2bNnY8mSJfjkk0/Qpk0bXH311Vi3bp3P+pMnT0ZMTIxWkpKSjHoKROHB7QaefRbo0wfxp44DAArtEXj06uF46vK7UOioOq2Hgu+XHi3w6n+vxZ5r4ypOvvoqcPXVQF5e6DpGdI7h8lvhRxIiOC/vAw88gK+++grff/89WrRoEdC1AwcOhCRJ+Pzzz73uczqdcDqd2nFRURGSkpJQWFjokWNLdE44dgy46y5g2TLt1I/JrfFExp34I6oRLCXyOctp+b+9tVjJ3zwt52+aS9RbeRcCLWdW23ZWlytbibqsFsxKrqxFPhbKUl1uq3Ks3LrVeuqxRVLOK+0p96vLWanH6qQq9WegUl5uNfmrUg0/7Sov5+X1WOqWuUo/3Tb51mXX3TZQlg9roOQfR8rR1ciG8udVbMOTAIDESPlbpKSIY2hhLUCnD47gshf3VyxPlpICfPklt7mleq+oqAgxMTEh+f2tPvaF7z8Gc4Pa37nQVezE3jtf4NjEAEGJyD700EP4/PPPsXr16oAHsQCQnp6Offv2+bzPbrcjOjraoxCdk/bsAS69VBvEuiQJ03r3w71D/4E/ohqFtm9UM5KEX+5sge8+aIPTzZWFaQ8dAi67DKgixYqIaoccPTVi1YJQP7P6y9DJXkIIPPTQQ1i6dCnWrFmD1NTUM2pn+/btSEhIqOXeEdUja9bIeZXHjwMA/oqMxJg77sKPSW0AAGYleqhGJNXtX03K9q+mcmXlAOXWpEZT1Rn7etrqBaj4c1hdtUCJzGq3Zj+rGmjHuralishr5WPf/VDr+DlWD/39EvFzvnJzkv41cHl2S1JfQy2SLN+61fTWUiXibJPvcJbJH7unSuWoz0mbfHvKakexRf75z4uj8c0nHXDVqN8Qs7MEKCwE+vUDFi2S32ciIgJg8ED2gQcewAcffIDPPvsMUVFRyM3NBQDExMQgIiICADBhwgQcOXIECxcuBABMnToVLVu2RIcOHVBaWor33nsPS5YswZIlS4zsKlH4+u9/5XQCZXOD3xIS8PcR9+Bo4yYwl4S4b3TGTsfZsGlRKi5+OBtx352U39/bbgMWLgTuuCPU3SOql4xa85XryBrH0IHszJkzAQBXXnmlx/n58+djxIgRAICcnBxkKYuzA0BpaSnGjx+PI0eOICIiAh06dMBXX32FAQMGGNlVovA0e7a8vJYSIf22U1s8fO+dKDZFABBauFNSoohavqmWd6qcN6uRWDUvVd2HVdlmVuiykNS1Xk2iIsKq5ZJ6bl8rvG7hcautt1rNdrPaOq2Vz5/p74bqvuZTI9fCx2OokVh9ZFaJwGqpwurzVV/jEvmHMrP84p80y8m1x6zyH/UNLKWIMMt/jNhNcmNmuxtrp7dBj8d/R9LS44DLBQwbJj/o0KE1eqpEVHP+9j+pjXbJGIanFlRnwYIFHsePPvooHn30UYN6RFSPvPEG8PDD2uHiy7vi8b8NgktYgRL5/54wq7e+B5HarVl/rA7ClPQA/X9ldRAqhNdEq4rH0rVh0Z/XHesGsD4HrvB/rvK1fvlLJajis0rAc8MGlTa01w+8da+t2ancoTxflzKQLVHem+NmZYKd5IZJlwNRprwxRc870Mv+O9otypVXpRg2DGjQALjxRr/9JiI6FwR1HVkiqiXTp3sMYmcP6I1/3XszXGoyLNUvJgk/TGyN3wYry3O5XMCQIcD69aHtF1E9wy1qw4+hEVkiMsA77wAPPqgdvnnzFZh6+9UwC/lLMZeQvKOE8DyuvAWr560aJdVFcpVtZKGmHGiTvSSvSVla22bPyV5q2/AbHfYzycvP9rNnRDcpzCsS6yMwq9XRHtezA/6iAd4RZpPHeZekTPqSHBXXKHeWupT0A2Utr8aWYgDA0ccbY9Cp7bjoyyNASQlwww3Apk1AmzZ+ekFEVL8xIksUTr76Crj3Xu1wxiB5EFvlzH6qP0wSPp90Mfb2bC4fHz8uD2aPHQtpt4jqDf0ugLVZyBCMyBKFi+3bgdtvl79WBvBu/+6YeudVypwqAdfZLFQo+YnEupRjixpGVOqLSvX1ObJqW2bdrb8IrD5yqZv05fNW91hnFaWtQuWUVW1lL2UCnD5R1ysqoOtTxYetZ81ypZ2TbqDcJd/nVCOy5fKEsAJrQwCAwyxPAnv52b546u9fIvn3Y8DevcCdd8qbJpgYmyCi8HLs2DHs27cPkiTh/PPPR+PGjQO6np96ROEgN1eOvBXLXzF/2b0Tnvtbf0Ziz1ElkTa89vK1KIxR0hK+/hp47rnQdoqoPjAqP5Y5sl4OHjyIAQMGIDY2Funp6ejevTtiY2MxYMAAHDp0qMbtMCJLVNeVlgK33gocPgwA2HF+Czz64E2AWYIEoeVVanmsbmjLY2m3QnerHmp5rWoeqO/IrNqMJKnLeVXaecDPNrFey23pVy3QR2bVeWq1kBsb6Ja0eh7X6VNkhb6S5HHeKzqg64P6oSspL4BLed/cbqBYiciWl8svRnGpvLvXcZscibWaXVo7h6Ib4/F/3YRpTy6C2S2AiROB3r2BPn1q8hSJiELmjz/+QK9evWA2m/H888+jXbt2EEJgz549ePPNN9GzZ0/89NNPiIuLq7YtDmSJ6rqxY4EffgAAHG0cg1H/dwdKbVaY4W/bLTpXbOnSErPv6o1/LFwnT8C7807g55+BZs1C3TWisCRvUWtMu1ThP//5D5o0aYLNmzdrG2SpHn74YXTr1g3PPfcc3njjjWrb4kCWqC77+GNgxgwAgNNqwahH7kRB4yhIgFckVrt1STC51LVPPddA1QKv+g9Vr1UMdFFSrb5n5LYyoYvEeq9O4HtVAm0N26pyYnV99Zcbq697tpFZ+IjMat3QvSbqsVvplMnvmrWe7anvkctlhkvZIthZJr9IZU75I7rYJkdizWbPP15cLhNeG3A1Om47it6/7gdycuQVLRYvrtnzIyIP3NkrOL788ktMmTLFaxALAA6HA5MmTcL48eM5kCUKa4cOeaxQ8O9hN+B/rZJgdsuDGm3gpHwdLcrlW6ncVDGAVXefcnve+uNvty1JHZSq530N0rQNAaoewHoPaNXbqlMKDP09UMVAt9pBsJZ3Id+YlMa8fnHp0zvU90R9j1wSJCWlwK0MZN1O+bbUKn9US2bPv0SE8gfLuDvvwPL/TEGTU8XARx/JqSi33VZNx4mIQiMnJwcXXXSR3/s7deqEI0eO1KgtTvYiqouEkAexRUUAgM+6XozFV3QNcaeorvozOhoTb7up4sSDDwJ//RW6DhGFK3VilhGFNM2aNUN5ebnf+8vKymqUHwswIktUN82eDXz7LQA5L/bJO24GYAJERRRO3TrVrURiRZkakZUgKZ8PJuVW0pbN0t2qh15LYHlOYHKrX59XtUqCPnKqTRzzbNN7W1xdtDcYEdgzUV2qgKSrp4StTeoWt1oqAXwfuyreJ3eZ8t5a5RdJWNTJd7oHV641OyV8feElWHXhz7hm704gLw949FHg7bcDeIJERMGRlpaGFStWoG3btj7v/+abb9C5c+catcWILFFdc/gw8H//px0+dtdtOOEjj4jIgyThub4344RdWZJr3jzgxx9D2yeiMKNO9jKiUIVx48bhrbfeQmFhodd9RUVFmDNnDsaMGVOjthiRJaprHnkEOHECAPDfrpdiXds2cuTNpS6FpYYylRtlgpCkRPFM5RIk9ZwuN7bayV663NiKKKQuf9XXp7Jui1q/y3D52fCg2mW2tPo1D9We8SSvM2jba1Uu7Vg+4VaeqKmKiKwaQXdblVvlE1pdBs3rtVHaMDvl278sjTAzvS8eXfuZ/B498ACweTM3SiCiOqV379749ddffd4XHR2NXbt21bgtDmSJ6pJVq4BPPgEA5EVF4fnrbgxxhyjcfHjxZRj0v8248FgOsG0b8MEHwF13hbpbROHBqO1kGZH1adWqVdi1axckSULbtm1x7bXXBtwGB7JEdUV5OfDww9rhK32vw0lHBCSXHLITZRXLawGVInpq9FWLzPrIjdVvjKCjXzpL3RChIlwqPG48NgHwtwSWV46s56331ra6Tp1Brmy1EdggRmi9VjNQIrPa6yB0UXMzINQVDNQcZy2P2LMt7TGV+srOtTCfBgAz3rjoBryx9i355JNPyisY2O1n8KyIiGrf4cOHceONN+KXX35BcnIyhBDIzs5Gx44d8dlnnyEpKanGbfH7JqK6Ys4cYPduAMCOxJZY2oWrFNCZ2RZ3ATbGt5EPDh0C3nwztB0iChNGbE9r1Nq04eyhhx5CVFQUMjMzsX//fvz+++84ePAgoqKi8OCDDwbUFiOyRHVBSQnw7LPa4Ut9bgCESZ7JXq5GLpW/O3XRUpMaxVNzZMukSmuTKre6vMyKhE5dP/Q5sl53VP9UvHNkdY+lz6X1d70fkhAVebL+nkco6HNj9ed1mzMINUKr5chK2qoEJj/5xPqVEdT33qREZC2n5TssJW7MadMf3XP3yuvaTp4M3H8/0LDhGT45onMI0wAM9+2332Lt2rU477zztHOJiYl4/fXX0bt374DaYkSWqC6YPl3elQnAt+d3xP8SWmrLMZnKJHlwqhRTmUkp+vNKKYdWJLdn8Zv/JUHZMUspJkkp8mBKmJViqkEx64p2n9KmuqSi7jHVY7+qyl3zd5/QlbNRXf+qeWxJLer74fIspnIBU6laAFOpPImrcrGc1hehK25YTrthLnEj0x6P1U06yn0pKNB2iCMiCjWTyYSysjKv82VlZTAFODmVA1miUCsuBl56CYC8ven0Xv1C3CGqLz5MuFxL1cVrr8mRfyLyq66kFqxbtw4DBw5EYmIiJEnCp59+Wu01a9euRVpaGhwOB1q1aoVZs2Z53L9gwQJIkuRVSnSfCzNmzEBqaiocDgfS0tKwfv36gPpeE/3798cDDzzgsTrB7t278Y9//AP9+gX2O5ADWaJQmztXXsAewKrUztjXJFGJxCqlVC2SXJy6op5X6usjfZKrIhroLzLpFRXVigRIPiKzVRTvSKvSRi1ubiMJoS1r5flEUKMIrN/rffDX74Cfj58IrcklF3X5LVM5YC4TcimVi8UpF69jpaiRWPW8qdQNU6kbhy1N8ENMO/nxc3Plf2tEVOedOnUKnTt3xps1zG/PzMzEgAED0Lt3b2zfvh2PP/44Hn74YSxZssSjXnR0NHJycjyKw+HQ7l+8eDHGjBmDJ554Atu3b0fv3r3Rv39/ZGVl1erzmzZtGqKjo9GxY0c0btwYTZo0QYcOHRAdHY033ngjoLaYI0sUSm43MG2adjj/oqtD2Bmqjz6Muwy9C+VJhHjlFWDUKMCs3yKMiAAYvvxWkbLtuMput8PuY0WR/v37o3///jVuftasWUhOTsbUqVMBAO3atcPWrVsxZcoU3HLLLVo9SZIQHx/vt51XX30V9957L+677z4AwNSpU7F8+XLMnDkTkydPrnF/qtOsWTN8++232L59O3bt2gUhBNq1a4e0tLSA22JEliiUvvoK2L8fALA54QL83ihRi9BpEVmlmEvlouVPqsdqJLZcLqbySpFYXW6sFpnV88qR1RfJd/G1pbhyn99IbIC5sV59FpXvEx7Ffxs1q+fzdfF1LpDzXn3xzJXV3pvK+cy6CK0WqVUismanrigRXDUSayqTC9xuHLA3xzabspRNZibw2Wc1f/5EVKuSkpIQExOjldoaHG7cuBEZGRke5/r27YutW7d65KKePHkSKSkpaNGiBa6//nps375du6+0tBTbtm3zaicjIwMbNmyolX7qdenSBXfeeSfuuuuuMxrEAozIEoWW8tczACxqd3no+kH12tLIzkgrzZYP3noLuPnm0HaIqM4KdFZnIO0C2dnZiI6O1s76isaeidzcXMTFxXmci4uLQ3l5OfLz85GQkIC2bdtiwYIF6NSpE4qKivD666+jV69e+Pnnn3HBBRcgPz8fLpfLZzu5ubln3cdnnnmmxnWffvrpGtflQJYoVH77DfjuOwBAdmQsNiS2heSuWErJpP7vVD5T3boF9jU+tjrVb4Qg6TZE0JaA0u98ql9lS78zrW75J9Rku1h/GyYEyKPPfpbdCijaCgT0FaL2Wuke289uvzV+jIol0YT3OfVYfT+1ZbeEcl65LZNvTWVu5bzyhiv3/+RIQY45GgmuImDFCmDvXuDCC6vuGBHVuujoaI+BbG2SdJ/H6hJ/6vn09HSkp6dr9/fq1QuXXHIJ3njjDUyrlOLmqx39uTPxme7boH379sHpdCI5ORkAkJWVBbvdjvPPPz+ggSxTC4hC5Z13tB8/T0mHkPjfkYwhJAlfRXasODFzZug6Q1SX6SeM1mYxUHx8vFfUNC8vDxaLBU2bNvV5jclkQrdu3bBv3z4AQGxsLMxms8929FHaM/HTTz9p5f7770fPnj2RlZWF/fv3Y//+/Th06BDS09Pxj3/8I6B2+ZuTKBSEAD74AABQLpmwMv5ieR3RcuGVG2sq1RX9/boc2UBWK9C6I0lK0a3xqs+V1eXQeq9y4KNoj+E7V9Z/p3z32Stf9kx+cVT7elTqa4DPI+D7K1fVvV/6XGeTS8DkEpWO1SKfh1spLqFFY+WGJaxo0B5OKJO8Fi4ESkv9vwBE56owHcj26NEDK1eu9Di3YsUKdO3aFVar1ec1Qgjs2LEDCQkJAACbzYa0tDSvdlauXImePXvWan+fffZZvPzyyx4TzxISEvDqq6/iueeeC6gtphYQhcL33wPKciY/NTkfx60NK1IKyuVPPHOp52hHnw7gL8VAndQFwGtHL58TvSq35WfDLD/f5J/ZZ3Og31D5eHD98/CbrnCWvzyE5J1K4K/p6l4bf/f7fU8qV1LTDrT3VfnBpTuvb0sNVZhMOGGKxA+WFFxVfgD46y9g2TJg0KAqHpyIQuXkyZPYr0wEBuTltXbs2IEmTZogOTkZEyZMwJEjR7Bw4UIAwKhRo/Dmm29i3LhxGDlyJDZu3Ii5c+fiww8/1Np45plnkJ6ejgsuuABFRUWYNm0aduzYgenTp2t1xo0bh2HDhqFr167o0aMHZs+ejaysLIwaNapWn9+xY8dQWFjodb6wsBAFBQUBtcWBLFEovP++9uPKxC4h7AidS761tpIHsgCweDEHskR6tbngtb7dAGzduhV9+vTRjseNGwcAGD58OBYsWICcnByPtV1TU1OxbNkyjB07FtOnT0diYiKmTZvmsfTW8ePH8fe//x25ubmIiYlBly5dsG7dOlx66aVancGDB6OgoACTJk1CTk4OOnbsiGXLliElJeVMn7lP1113HUaOHInXX38dPXv2hBACGzduxEMPPYTrrrsuoLYkIQKdHVG3FRUVISYmBoWFhYYlVBOdFZcLiI8H8vNRYrLg5sufRInFjtIo+e/K8gZyGK1cmczqtskfgG7lW2GhLgGqRk/1k5AArwlgXhPC1Lr6Yy2yV02kr/K31gZ9glT5uV/D3wlG/D7y+3z9nNe/tv7Oy0ty6V53ZZKXel59/0z6SV5KFN9UKl9gLimXj0/LYX7plLxzj/n4cXzw17uIFk6gQQPgjz+Ahg39PCGi4Arl72/1sZOmPwNThKP6CwLkPl2C7Aee5thEUVhYiAceeACLFy+GyyV/bplMJgwdOhRvvvkmYmJiatwWc2SJgm3TJiA/HwCwNfoClFhqZ/kVouq4JDO+t7eSD4qLgS+/DG2HiOoYIYwrVCEmJgbvvfcejhw5gnXr1mHt2rU4evQo3n333YAGsQBTC4iC7+uvtR83xVyoLaFUEVVTcmRNugxVi3wslKicUP4M9bkqir/on57fnFhJORa+mqsIiIqqJy6dDa+c4Mqqy/Wtpg9nE6n1+3z91ddV0J6X7jpRuU51jfq7X/nHIEzyPw5hkW8lq/JRb7dhLdpiQImy09dHHwFDhlTzIERExmjevDmaN29+Vm1wIEsUbMrasQCwLaZ1CDtC56KdtkQchwONUCKvKet0ArW0KDtR2DNqhQFGZD2Ul5dj3rx5+O677/Dnn3/CrS2ULlu9enWN2+JAliiYTpwAtmwBAGTZY/GXpWHFovblnpFZIac5QkiekVnhVqJufnJlq+IvX9Mf9bHVyCx0UVKPh6xmJYEzjdD6uq7GqxRUEy2uUWQ20HxcfTRcX0933uNYjahqOx8o97k9w7na629SI7DKefVYjd6Xy/9IhE3+qJciHHAD2OJIwbUle4BTp4D164FrrqnZkyQiqgVjxozBO++8g+uuuw4XXXTRWW24wIEsUTCtXQuUyyPUHVEtQ9sXOmdtsyfJA1lA/oaAA1kiWR1ZtaC+W7x4MRYvXowBAwacdVscyBIFU6WvS7Y3bAW4K7YZVWeiC/W2zDNmJwl19QLlfjWqaNKFRwP5vKwu11SLXOois5Xu9rfOqr8IbW2sclDj3Fd/fdO14/d3TA1eS+21Ud8Q/Xqz/hbj1UfH3fAfgTXp6qqRV/Ux1WN1ZQu3kiNrg2c95fbnRhcAhavkcwF8hUdU32mbkhjQLlWwWCxo3bp2Uuu4agFRMH3/vfbjLw3l/aXVKa1SuWdRd/pSd+zSdv5Sd3MqV4t6Xtn1qdJOUIHuNFP9rlRSpVQH+T6hK/6u9fsYtcjvL6FA89589L1i9zPPUt39Wj0JHsXnLmrqeW03NUkuZnmQqh3r77fIxa0Wq8mz2C1ycdjgdthwLKYpsqxN5I5v2QL4WJiciMgoY8eOxWuvveaVG3smGJElChanE9i+HQCQbW2CE5YGIe4Qnct2NEhBcuFf8rrGP/wA1MJXfERhj5O9gmLTpk347rvv8M0336Bjx45e2+guXbq0xm1xIEsULDt3AmXyAvV7HQkV57UF8ZUUApfnrZo54FZChCblArdZt3aW9rWzqJRm4GtW1hlERP2lGhi5OGJNmg40ZaCaVAOf3fA3CSHA11BrR8tA8FzaTEKlFAJ3pZOVb9UJftClFug6o2YmaMtwKf+IJHUSmMWE/zVujRsK5T+ssH49B7JEFDSNGjXCzTffXCttcSBLFCw//6z9uN8RF8KOEAG/RiZVHGzeHLqOENUlnOwVFPPmzau1tjiQJQqWnTu1Hw85mntnqKuROnXbWF1EFn4is0KdBFa5LfUadVkmXcSuVj9SvaK1yml/k5uqE0iQt5rltgJaZqumzjK6rb0cuggthNDee+E1e8H3i6j+W5CUJ6ouYSOpk8bM6iQwNcqvHFtN+NMai3xbFGJLTwDbtsm52mexBA4RUShwIEsULHv2aD9m25uGsCNEsv0NExH71x55stfhw0BSUvUXEdVnzJENilatWlVKjapaZmZmlfdzIEsULPv2AQBOm6z4M6IRhFlJePSKIirRM7cWgpWP1Qgt1NPK/Rb1vJo36SOiJ2mJmXIdr71Qz4yQJGPzZANV0xzYM8iVrTYSW9ONE/xcJiBpEXS/kVl9zrP6vqqRVjUSq3ZO/bekHKvtut0CktWEAzHxSP9L+QPr1185kCWioBgzZkyttcWBLFEwlJcDBw8CAHJtjfkVLtUJWQ2bVRzs2wf06xe6zhDVBYzIBsXDDz9ca21xIEsUDFlZ2ooFRxxNIWwWCIsyo9zsZ1CrRjpd8o2ki8apeZEo94zM+gwNKpG6ihq+c2fPKkJbXa6sTm1ukFCdM82V9btigUfjNbzfz9vjsWqBeqe/yKzXJhOeObHQRcclf09YSJBcwKGmzSvO/f57FU+C6BzBgWzQ/Pzzz/j999/RtWtXJCcnn3E7HMgSBcOBA9qPhxvGwhVhqbQjl/q1sG5Cln5ZLZdneoCaSuCG59fL8gV+BjD+BrT6iUdaJ/w+o8AFOumrjqp2MOzv/mqWAqv8rukHtB67fwFeEwUrxq++/zDR+lz5bgFkiyYVx1lZfjpORFS7Xn/9dYwbNw4WiwVmsxlfffUV+vTpg2nTpsHlcmHs2LE1bos7exEFg5JWAAB/2BuFrBtEleU3iAaUtWZx6FBoO0NUF+i3KqzNQpqXX34Zr732GpxOJ0aPHo3JkycDADp37oz58+cH1BYjskTBkJ2t/Zgb1RQuh8V/hE5PeH4vLmnrbHkuuaT/klq+RLcIvy6ip4/Mai37i9AGQz2J3FbLx9vmvWKZ8p6r76f+ffNKNfBsQGjRfuXYXHG/kIAyWJAbEY34U8c9/o0SERnp+PHjGDhwIADg9ttvx6JFiwAALVu2xIFK32DWBCOyRMFw9Kj245+O6BB2hMjTn5HKv8c//9TyuInOVZIwrlCFyy+/HN9//z0AoEmTJigqKgIgL7XVpEmTqi71YvhAdsaMGUhNTYXD4UBaWhrWr19fZf21a9ciLS0NDocDrVq1wqxZs4zuIpHx8vK0H/MbxaA80gyXwwSXwwS3TSkWCW6LBGGSI2j6b6O0D0O3kIsyKUFyC6WgoriUoruv4hqlqG3qJjhIQgRtWS2/37pJqOWdG/wL6S8aH89Te02UIiRJLia5aOeVfytuswS3WYJQivpvyW0D3DbAZfcs5Q7ApZT8qKiKB/7zz+A9byI6Z9155514/PHH8dxzz2H9+vUoLy/Hxx9/jNGjR2uR2poydCC7ePFijBkzBk888QS2b9+O3r17o3///sjyM6kgMzMTAwYMQO/evbF9+3Y8/vjjePjhh7FkyRIju0lkvPx87cfjtsgQdoTI03FHpX+Px46FriNEdYEwsJDm7rvvxtGjR/H000/jvvvu03Jlr7rqKrz88ssBtWVojuyrr76Ke++9F/fddx8AYOrUqVi+fDlmzpypJfZWNmvWLCQnJ2Pq1KkAgHbt2mHr1q2YMmUKbrnlFiO7SmSsv/4CAJy02nE62goAMJcqGxmUq7dKXeVWy2fVR0bVfFd1FQN9sqVUcU6fT6ue1xbO95MzC6322X/6nvWWtSHMmZWE8FqCy5Btb7XGldtqVo/Qcl/1ObEm3a3yCe9WcmPdlorzar5sUaSjouHCwjPtORFRjR3T/dFstVoRERFxRm0ZFpEtLS3Ftm3bkJGR4XE+IyMDGzZs8HnNxo0bver37dsXW7duRZmf3C2n04mioiKPQlTnKP8uT1kd1VQkCq7TVlvFQXFx6DpCROeM6Ohoj3Kmg1jAwIhsfn4+XC4X4uLiPM7HxcUhNzfX5zW5ubk+65eXlyM/Px8JCQle10yePBnPPPNM7XWcyAgnTwIATtntKGuorP/qlO8y627V5WRNZeqi+H4igrpjz7VGlWtd+rBhNZFZfYQPuvMGqnaDBN0aqH7vq6Lts+JnQwMfSw148vfYvp6nv1UI/NFHYs2et1ok1qrc2iqO1TolEZV+DZSUVPOARPWbBGPy5bn4lnEMn+wl6b6WE0J4nauuvq/zqgkTJqCwsFAr2VxChuqi06cBACV2K0obAqUNgbJICWWREsodcnHZ5aJN+lKLqWKgUpk2QcklAFflyV76n0WlyV+ek7y8JoHpJn95TQY7i1wvr0ldfiZz1WjpRUlXalsNnqfXBDF/r1FN8+UCeV11k730xa0Wi1KsnsVlA1wOAZdDoNxW6QV0u/0+JBFRXWRYRDY2NhZms9kr+pqXl+cVdVXFx8f7rG+xWNC0aVOf19jtdtjt9trpNJER3G5tWaNSM5dupjosSCtVENVZRm1ewA0RDGNYRNZmsyEtLQ0rV670OL9y5Ur07NnT5zU9evTwqr9ixQp07doVVqvVqK4SGatSlKvcKqEsGnJpIJfyBhLKG1SKyNqUoltSyevzVQhACK9lueQIq3KfPjKrPw40Mqs9NvxHFgNVTXT1bDbHqY2lvbyWIvNaqszPWpEBzmiu0fqT+tdKH5k164pFLvqIrNsu4FKKyeSqaJ+fs0QUZgxNLRg3bhzefvttzJs3D7t378bYsWORlZWFUaNGAZDTAu6++26t/qhRo3Do0CGMGzcOu3fvxrx58zB37lyMHz/eyG4SGavSIEg/A54o1Oxl5ZUO+O0WneNqa6mtWkrJouoZ+j3n4MGDUVBQgEmTJiEnJwcdO3bEsmXLkJKSAgDIycnxWFM2NTUVy5Ytw9ixYzF9+nQkJiZi2rRpXHqLwpvZrP1oMrlQFuO55pUacVMnZqkTryQlUFZxv26Wl0rNI1cmhQkIP9vWwkcUUrlfnWGmn0jmuxVdZFjfZtX0EVKviRXVbN0b1G/oBCpeOt3X7sJPSLfG3avJ66abACd0572W31Kb1i+3pU36kmsIiwBs8j+0BuWVJng1bFjDzhMRnbnU1NRq6wghcPDgwWrrGZ6wN3r0aIwePdrnfQsWLPA6d8UVV+Cnn34yuFdEQSRJgMkEuN2wuDiZhuqWmOLTFQeNGoWsH0R1glHRU0ZkPWRlZWHSpEmIqryzYCUnT57Ek08+WaO2OPOEyGiSBDgcQHExIl1OuBvJX+WWCTkfUY2kSso3vKZySblVLlcisZJZDbspx1VGZv1EYt362Kp2sed5k8+zVUdmz1DAEdrKHQlUNf31udmBn/743cK3Fn9hCb+RdB31ft2GCOqxFplVI7I2txaRbXLqVEU7fibVEhHVtvvuu8/v5P+8vDwOZInqlOhoeSB72hnqnhB5iDuubCJjtQJNmoS2M0Qh5jXBshbbpQo2m83vRleAvKmWw1GzDYQ4kCUKhkaNgNxcxJw6jcgY+avck0rktbxU/m9oKpWPXcr/7YqIrBKhVXNm1VxZyXeEVkBoIUXtw1MXiRX60KOWXOl53m8kttJxlWu9VsVfQNPPdV6z9w1U+bG8+nOWv5Bq8gtN6KK+wuT7/dDXF34is8IsPI8tAiaLHJFNKFC2pU1MlFNgiM5lTC0Iivj4eGRmZqJFixY+7z948KDfaK0eP7WIgqFZMwBAg9NlaOkoQGJMIWxRpbBFlaI80i2XCKA8AnDZlWKT4LJVbJDgNss7NfldjkvhsWyTbrMEv8ty6TdMUM77XX7Lx+MFzM9SUjVdhutMl+QKuJvVLI2lX8LsbIr2enttXKG8H9qLAd+/GNXluNSibqahFGEWEGYBySwPZKNLi9HolJIj27KlsS8kEZEiPT0d7777rt/73333XVx66aU1aosRWaJgqPSXZZOCU8hNjAlhZ4hk5/15vOJAWU2G6JzGiGxQ/OMf/8BVV12F+Ph4jBkzBk2UtKZjx47h9ddfx9tvv41Vq1bVqC0OZImCITlZ+/HioiPYc6ELhSURAIC80/KkL1eJSblVUglK5fpm5datTgLTpxj4oqYbqKkESjTPa1kuLT1BOa0uqqCrpj2S8HG3URHRqtoNcEmugCPGNahfbZtVRLD9N+rnWP3uzO2ZalBtO+rbq5v8BUnAZBJo/cefFddceGHVbRIR1ZLLL78cL7zwAh5//HH85z//QbNmzSBJEvLy8mAymfDcc8+hT58+NWqLA1miYGjVSvuxWXYR9nSND2FniGTnZ+dVHLRpE7qOENURnOwVPOPHj8ctt9yCTz/9FAcOHIAQAi1btsSNN96ICy64oMbtcCBLFAxt22o/djhyFAWRUfijRF4/7/gpOTJbekpetb7cId+alXXqXVYlEqtMAnOb9ctxKQ2rUddKH8RCC6nqJn/pIrDVTv6SPK8/k7zU2trVTBLCf7S2mglkZxJF9V7mzF+/qq5Xk19k+mXNtEiq/n1T3xY/kWmvyV/6TkiAJAm0O5Rbcd9FF1XfQSKiWpSamoqxY8eeVRscyBIFQ8eO2o/NfjsRwo4QVbho/xH5h6gooHXr0HaGqC4wagZpULckDA85OTmYPn06fvrpJ5hMJlxyySUYPXo04uMD+8aSA1miYIiLA+LjgdxcxO8pwvm2HOxr0BwAcLRBNADgzwY2AICrWMmVPS1/8JmVpWfdSq6sUCKx6raj/pbjqkyLpLp915X0kdgAcmWr2yTBbyQ20M91LRrsfaG2OUENt7cNKHpaXWTVT3RU/z4E9NWi5PWDx6F+y1rJXx/8EUD8n4WI/0tZQ7ZbN4+tlImIjPT777+jV69eaNKkCTp06IDPP/8cLpcLF110EdasWYP27dvXuC0uv0UULF26AABsx12IOOx/IWiiYOj5y4GKg8svD11HiOoSYWAhzb/+9S9cfvnl2LlzJ1566SXYbDZ8/fXXuP/++/F///d/AbXFiCxRsHTrBnz9NQCg7c9/IOmavwAAvzeIBQAci2gAAHDZ5f+WLoeyQYKSK2uyKRFaZQysbV1rUcJ0ZX6ikpXpI3dqZE+fE2tgruwZr3JQxVa1+iitV4S2hr9EfEVbvfKK9XX9bhms1HP7Pu/94JV+1kIMnq+70Pelmnxdrz65K5JnL9uxv+L81Vf76RQRUe1bvXo1vvnmG5hMJohK317dfffdmDp1akBtcSBLFCyVol7NNpYisf9xAECsQ97rPidCTjEojFBSDE7Joxm3fAi3vEoXXMqtuuOXOlAyqxO34D3Iqo7XYK26FINaGMCeTcqY1l8/qQRCG3BX8/z93V15IOg1+NcNXPUDXLfn+YDeC33ahHqotCmZ1Db9pIT4mxSm9km5NZW6cf2uLPmgcWMgPb36vhGdA7hqQXA4nU7ExsZ6nS8uLkZkZGRAbTG1gChYevYE7HYAgHVNSaVtY4mCq+v+Q0ChsjVtRgZgYUyDCABTC4IkJSUF+/bt8zh3+PBh/Otf/0JGRkZAbXEgSxQsERHANdcAAEy5bpy/6w80sxShibUYTazFiLSXItJeCpPdBZPdBbdDwO0Q2pa1bptSrJJclK1r1S1r1S1sUXn7WkkCpMrHqN0NDM7gw9lrUnB1W9X62LrW7xa1/p5fNc9bi8IopfJxxc8CEMK7rrqdrJ9tfisKalB0dZU21ceu7hek321zlXZM5RJM5RIytu2quGjgQP8vDBGRAfr164fFixdrx8XFxUhOTkZZWRlee+21gNrin+FEwXTTTcBXXwEAzvusEIUXNwhxh+hcY3K7cf3Wn+UDiwXo3z+0HSKqSwxKLWBE1tMrr7yi5cYmJCRg2bJlaN26Nc4///yA22JEliiYbr1VSy9I/KoQkS4noi2nEW05jQbWUjSwlsJid8Fid8FtcytFjsS6dMVtlYt2rERo3WZAWCS5mCAvqq9GJE0SYJIMWyqxMiFJcq7qGUZJq72u0n3+no/Wh0BVisbqI7Fa5DXgCKz+fl3xiLJ6RmC9Iq36X7b6CK1bLtpjl0tyKZNw2a59iC9Ult0aMABQ9jgnIgomdcv0iIgI9O3b94wGsQAjskTBFRMjDx6WLoUt34Um3xcDl4S6U3QuuXPDhoqD4cND1xGiusiofFZGZD1cddVVHqsV6K1evbrGbXEgSxRsI0YAS5cCAOLfK4S1qwsA4DCXAwAsFvm41CZPMXfZ5C9OzMryW+oqBq5yz1ULKi/zZHIpP/hbfsrfFqa1wU+bfh+ruj5UsexWdc+vRqsS+DiuvNKAds6lq1tpS2DfbVa9uoF3X4UcMa9cV3ge1vh5qP8m5H9SWrPxxwvRZ/du+SAxEbjhBn8tExEZ5uKLL/Y4Lisrw//+9z/873//w9133x1QWxzIEgVb//5AixbA4cNo9u0pNM08iYLUhqHuFZ0Dbt+yCRa3Msq95x6uVkCkx4hsULz66qs+zz/77LM4efJkQG0xR5Yo2KxW4JFHtMOec/ejzG2GxeSGxeSG1eKC1eKCZBWQrALCJhd9TqxWrJJcbHJx2ypWNNCvbKAVk77IubR+VwOoTbW9ckJN6Fcg0HdJtxKBR16stnKA0kalnFbP1QuqyZH1l+vqrhRN169OcIbPr2KVArmYnUDDE07c+eMPcn2TCfj738/stSQiw61btw4DBw5EYmIiJEnCp59+Wu01a9euRVpaGhwOB1q1aoVZs2Z53D9nzhz07t0bjRs3RuPGjXHNNddg8+bNHnUmTpwISZI8Snx8fG0+tSrdcccdePvttwO6hgNZolD4+9/lhegBdP7qMMQRoNRtRqnbDEkScjG5IZnccFsE3BZRMalLV9TludSBrDyo1U0AqzQRrHLRBrC6ga062Kx2gFuDCVd++RnQ6h/LZ7vVDYa1gZ2oelOEapavqjyByt99fge26rFbVF0q90X3/Gr6muoHxSaXUpxyMZ8GBm/ZhCbF8uYbuP12ICmp+oaJzjF+/+CshRKIU6dOoXPnznjzzTdrVD8zMxMDBgxA7969sX37djz++ON4+OGHsWTJEq3OmjVrMHToUKxevRobN25EcnIyMjIycOTIEY+2OnTogJycHK388ssvgXX+LGzYsAE2my2ga/i9ElEoREcD//gH8PzzsJS7MWD2/7B9LAcWZAxbeRn+tuW7ihNPPBG6zhBRtfr374/+ASyNN2vWLCQnJ2vbu7Zr1w5bt27FlClTcMsttwAA3n//fY9r5syZg48//hjffvutR16qxWIxPAp70003eRwLIZCTk4OtW7fi3//+d0BtcSBLFCrjxgFvvgkUFaHH5wcQe90JZLaMRbnLLN+vRuIs8p/ybpt86ypTJnmpk48qTfKS75BTBQBAqE251IlHnpOJJN216nVek8HUJax8TRKrpTQBo5cDkx9EvvGe1KXees6uUlMC5J/191Uz2UtX34vPKLN8suJ1h8ex/tbrebkq+l35/O07NqFZ8Qn54OabgY4d/XSKiIxUVFTkcWy322FXlmQ8Gxs3bvTaEatv376YO3cuysrKYLVava4pLi5GWVkZmuiW4Nu3bx8SExNht9vRvXt3PP/882jVqtVZ97Gyxso3kiqTyYT27dvj+eefx9VXXx1QWxzIEoVK06bAY48Bjz8Os1vggdlrMP75W0PdK6pnGjpP4+/bVlScYDSWyD+DJ3sl6VJ6nn76aUycOPGsm8/NzUVcXJzHubi4OJSXlyM/Px8JCQle1zz22GM477zzcI2y4yQAdO/eHQsXLsSFF16IP/74A8899xx69uyJX3/9FU2bNj3rfqrmzZtXa21xIEsUSmPHArNmAVlZSN+SiY5rDuPbzu0BAMKti4JqkVn5hEsXddPqmQCzWVmqS11+SbnVlmXyF0VUT/uKFsJ/pNDzXBXR2yqofdDX95lbVtNfNH6en/56fYTWI6Kpj9L6ieZ6RWL9UZ+fyTP6KqSKiLhQ7zPrjnURc61Jt+9bU7nAyC2r0LhEyY0dMgS4hAsXE4VKdnY2oqOjtePaiMaqJN3mL+o6rfrzAPDSSy/hww8/xJo1a+BwOLTzldMZOnXqhB49eqB169Z45513MG7cuFrr62effYbjx49juLKW9eHDh/Hf//4XSUlJuPXWwAI6nOxFFEoOBzB5snb41OxlaFDiDGGHqD5JOPEXhu5aBwBwmizA88+HuEdEdZvRk72io6M9Sm0NZOPj45Gbm+txLi8vDxaLxSuSOmXKFDz//PNYsWIFLrrooirbjYyMRKdOnbBv375a6afqhRdegMsl50I5nU706tULs2fPxn333RdwhJoDWaJQGzoUUHKCzvuzEA+/txplp6wQZSaIMlPFV13KTHa3VchFXbXAoSt2oNwhl4rzklzsStFvc6sUt0UuwqwrpipKLS/VVe1MX6H72ddKAvrVCnT1/C23pd/q1fd2sKIiSlujJ6Qr6jbB+tex8rJoFkBYALdZgtsseb8Pyta7XqsVlAuYygXMTgFziRuP/rAUNrf8y8L+f/8EUlNr3m8iChs9evTAypUrPc6tWLECXbt29ciPffnll/Hss8/im2++QdeuXatt1+l0Yvfu3T5TE87Gb7/9hu7duwMAVq5cCSEEdu7ciY8++gjz588PqC0OZIlCTZKAWbNQYpUzfUas3ICLDmSHuFMU7vpk/YLLcnbJB3FxwIQJoe0QUbgQBpQAnTx5Ejt27MCOHTsAyMtr7dixA1lZWQCACRMmeKw0MGrUKBw6dAjjxo3D7t27MW/ePMydOxfjx4/X6rz00kt48sknMW/ePLRs2RK5ubnIzc312IBg/PjxWLt2LTIzM/Hjjz/i1ltvRVFRkZYCUFtcLhcaNGgAAFi1ahX69esHs9mM9u3b448//gioLQ5kieqC88+HY9KzAACzEJg66yNE/lUO0ykzUC4B5VLFmqoWAWERcNvkoq0jq0RfyxsA5RFKUSKz+mMtMqtFaCutQWuVtMisvugjg/Jas0oUUYkS1njDg+o+8Gvwy8FvJFY91m1G4HW9bg1Y7/rC+7Gqo0Rc9UXdjKJiDV/5datY41fSIuT6SKz+Gv1rLLkEJJeAuRQwlwKNj53C//24tKLCG28AMTE16DwR1QVbt25Fly5d0KVLFwDAuHHj0KVLF21pqpycHG1QCwCpqalYtmwZ1qxZg4svvhjPPvsspk2bpi29BQAzZsxAaWkpbr31ViQkJGhlypQpWp3Dhw9j6NChaNOmDW6++WbYbDZs2rQJKSkptfr8OnXqhLlz52Lv3r3473//iwEDBgCQJ63FxsYG1JYkRCDfkdV9RUVFiImJQWFhoUdCNVGdV1YG9OgBbNsGAHj/0l6YdP0tKG8g/xdVl9/Sk5RJYZI6satMgqkMnue0+5TzumWavI51t/4mR1WlxukGfpaS8lu90oBX8jMRy++xbrKb15JaLh/3e7UlfLbt3VHlbt0SZvoJXW5lyq06yJV/lm/dZj+TvHTPx1wqn7Cekk88+uPHuPHwj/Kd110HfPGF9vhEdVUof3+rj33+v56H2e6o/oIAuZwl2P/i4xybKNasWYPrr78excXF6NGjB9atWwez2Yy5c+fil19+0dbDrQmuWkBUV1itwPvvA126AKdP487NP+D789tgxSUdQt0zCiPdc3+rGMQ2aCCvVcxBLFGNnMkuXDVtlypceeWVyM7OxqFDh9CpUyeYzfJf7/fee2/AbXEgS1SXtGkDvPIKMHo0AGDykkXY3eyfyI1pjPKGchWXQ/5EFGbPW1iU6J1VwF2uRmnlWzUSqy3DpV+WS71VIpImNRKrj9CqAczKUcrqIpPV0G/K4MUrGiz8R4j9RF69ls6qJhJbk8lc+qXIKjrjZ1MD3RJalSOxgJI+YNFFbc2+H0N7n5SNLtSIbOyxE5iw9b8VFV9+GWjZstrnQkQUbI0bN/baGOFMcCBLVNeMGgWsWAF8+ikaOYvxyifv4e7ho1HOlHaqihB49Jf/oplT2TkoI0PeBpmIau4MJ2fVqF3StGrVClVltmZmZqKgoABdu3ZFZmZmlW1xIEtU10gSMG8esH07cOgQ0rIz8dgXX+I/g24EUBGlcykbJAirZxjSJSRILiWEp0Zk1QitGpktU88rl+oitNpGCtXk0grhP/pZbV6tn88wSf/h5qudmubEVpPX6i8SW2X+q37RcX/byXptaqDcVorAApWisObKubGebWoPreuv+n5Zi9247ff16PXnbvlEbCwwfz5TCoioThozZky1dSIjIzF27Nhq63EgS1QXNW4MLFoE9O4NlJdj2M9rsaN1Er7qxF2ZyFvbY9kYtWtZxYl33gESE0PXIaIwxRzZ4Hj44YerreNwOGpUjwNZoroqPR2YNk3Ll33u88U4ammGbZ3kvbrVrWrhUEKiFvlWMlXkkAq3HAZ0aRFZ+dhdqubOqud959IK/eoHSlTRI0KrBv30Kx2o/ERi9SsOqPxFQyuf9xfdrTYC6666vhepoq7QZXZ4rUagRmK1lQd858RqKxFox9COtaitn216ta1n1dzYEoHGJSfw3OaFsAolXPvoo4CylA0RUTjJz89Ht27dqk0nqIxJd0R12ahRwN/+BgCIKC/D61/NRfOiwhB3iuoKi7scz218F82dyr+Jnj2B554LbaeIwpmP9aZDtSlCffbll1+iTZs2sNvtMJlMWmnevDmysrK045pgRJaoLpMkYOZMYM8eYMMGND9VhFkL5+Fvgx/AXw3kbQfL5c1RYLHK4TqLrRwmJYTnVkJ7bpcSmVUjsmVyGNBVKh9LTvW8EplVIrbKiihaRFboVz+AvN6/z667dCd0KwbUPL9V165AtbmsfiOvWl88T1T1tV/FGq5+ViPQcmA9c1+91oJVI7C6W4+8WP3ntpoLW+Z5bC6Vn8OETf9Fl/wD8snERGDJEnkZNyKiOuyf//wnrr32Wlx99dXa0lsAUFhYiOHDh+PTTz+tcVscyBLVdXY78MknwKWXAllZ6PhHNl789AM8OOJuuE0muCKVwaiyZpbDVga7RR5FqgNal7JpQqlL/sBwlsqDnbJS+SPAZZfPixK5LWFRbz0HttpyUKXyjSQBUNMO1P6KSveh8mQwz8Gn12DTa9JYFekBZztA1Z/XDcYrf7WvDlBrPHC1eG50oA1YfaQSeFxnrngM7bUo9+yvOqC1nnJjxO6V6H9E3jwDDof8byQ+3ufTJaIa4qoFQXHw4EE89dRTiIuL8zifl5cHALjhhhtq3BZTC4jCQVwc8OWXQEN5Mdmrsv6H/1v5RYg7RaFyfeaPuG/3iooTCxcC3buHrkNE9YS2XbUBhSqcd955sNvtXufNZjNaBrj2NSOyROGiUyfgo4+AgQMBlwt/27gWBaYYzLjlCgCAO1r+u9RidiPKXgIAcJjLPZooVcKBp8vliGxxqU2+LZFvS63yR4LbpqQenFYitWYlQquLTnqkFahRQ3XCktujqnfk1etWicC6dPV1KQm+7qs4X8PfFvoltLTznvcLyXvyVnWpAy7lm31/qQRuq+f1HtvP6l4js24ZNGuxwJVZ/8P/bV9S0fmXXwZuu63Kp0tEVJccOHDA5/mmTZti9erVAbXFiCxROOnfX86ZVYz/4XPctGVrCDtEwdT96B489/17MKsj3nHjgPHjQ9spovqEk70M87e//Q2nTp3yeV9ZWRn++9//ol+/fmjVqlVA7TIiSxRuRo4Ejh4FJk4EALy4aDHKyuz4fEA7AIAQEhpY5GTKJvZiAIBdmZ1VpoT/SpSwYVGpAwBQaJNvT5TIX/WctioRWpNcz60m42vrb1X8DaxtW6vb1tZrMX99RFZZQkqbOKaPzGo5tN65sl6RWLef3xK6mWgVE7Xk+sLP/dqtpVLEVcuBVW+VCK0SYdVudZFXfQRWP8lLmCqen7qRhcmpdF/JibUUA92O7MVLa+bD6lbCsyNGyNFYIqIwsGjRIowePRrdunXTzv322294++238e6778LlcmHo0KF4/vnnA2qXA1micPTvfwP5+cCbb8Ii3Jjyybs4mTgM313cLtQ9IwNcemQvXlvxNhzqlms33wzMmVPxhwUR1Q5O9jLMTTfdhBtuuAEPPvggEhMTMXfuXGzatAnXXnst3njjDdx0002wnsGqKxzIEoUjSQJefx04eRJYsAA2twszp7+L0X8bgR/6JwNRcrVIixzaa2yRI7NWJdmyRAkTFimR2GOl8hpef1mVW7N8e0KJFjrVyKyk+8gQJq/ta9WAodfuqLqIrBaJVSOz2va3uoisSx+RFdVOnNAir27PY0lbeUD9QU3sVTcv0C+VJdU4Auuywff9Ws6sEgVWV34wqaFspS+uShtPqJHY00CfA//Di8vfhU19YW+8EfjwQ8DCj28iCh/vv/8+5s+fj9dffx07d+6E1WrFK6+8gtGjR5/RAFbFP+eJwpXJBLz9NjB0KADA7nJh5vwFuObH3SHuGNWWG3dvxpRv3qkYxN50kzzhz2YLbceI6imuWmAcSZJwzz334Oeff8aWLVtw33334ZlnnsF5552HRx55BD///PMZtcs/6YnCmdksL70kBLBoEWwuF6ZO+QgvSn2x8tr2QKRcLcZyGgDQ0FzicXmxkgubb5FDuJEWeYFYmxIeNStRw0LlU1i92q2FGwXc6va3yqeJtni/3xxZz4irqVwXkS0XunrqagZqBLOK10PLcdWte6uupKDfvcHP9rKV14J12TwjseqqBG6b7lYfibUpEVglEquudqD/jaZGsqUySd7oAID1hMA/Ni3HqM2VltgaNgyYO5cbHhBR2LvkkktwySWX4JVXXsHSpUsxd+5cpKWloXPnzhgxYgQeeuihGrdlWET24MGDuPfee5GamoqIiAi0bt0aTz/9NEpLS6u8bsSIEZAkyaOkp6cb1U2i8GexAO++C9x9t3zoduOJF79G/3d/wYkyB06UObSqUabTiDKdRrzlOOItx5Fk/QtJ1r+Qav8TqfY/kez4C8mOv3Beg+M4r8FxxEWeQFzkCcQ0KEFMgxLYGpTB1qAMiHABES64bW64bQJum4CwqJOj5CJMSpGUgaJuBq/kEpBcAiaXPIiVytyQytwwleqK0yUX7ZzLu5S5YSpzQypXiksupnK5SG554KyPjKh9E2YJwiwpqQRyOoE6iHXZ5LQBl10pDqVEKEU5Lm8gUN5AwNXALReHXNx2AbddADa3XCxCLgqpXIJULsFyWoK1EGhYUIqXPnvPcxD7yCPAggUcxBIZjasWBJXdbseQIUOwcuVKHDhwADfeeCNeeeWVgNowLCL722+/we1246233sL555+PnTt3YuTIkTh16hSmTJlS5bX9+vXD/PnztWMbv0YjqprFAsyfD0RFAdOnAwD+9c4KLD/dHu8/woXyw0XciWN4/fP56JB3WD4hScCrrwJjxoS0X0TnCqPSAJhaUL3k5GT8+9//xr///e+ArjNsINuvXz/069dPO27VqhX27NmDmTNnVjuQtdvtiOdWi0SBMZmAN94AmjcHnn4aAND3o12IOliCr1/tBGdDK+KshQAAh/L9f5SSLKAem7UZVjK3slxXuVu+LVNuXeXKcbmkpRa4lW1sTeryUn6+71G/Sjepk7jU1IIydTkut3JeuVWX1tImfVX8RlBTCLTvltRZ/Bbd/Vp9z75pmxxYdekDdvlYjcYCctQVqEglcNnVYyV1wKr0y+pW2lb6q98gQnntTGXKBhan5Me67Od9mPLlQjR2KussRkbKkfabbgIREfkW1MlehYWFaNKkSbX11qxZg+bNm+PCCy/EyJEjtb13fXE6nSgqKvIoROcsSZKX5po7V5vV3nPzATw47Ds0y+T/jbrI5HZj1NpVmP3JrIpBbKtWwMaNHMQSBRtTC8JO0CZ7/f7773jjjTeqzX3o378/brvtNqSkpCAzMxNPPfUUrrrqKmzbts3nvryTJ0/GM888Y1S3icLTPfcAqanALbcAx44hLvMEHrhzNVZM6oDd/RMQZZInfznMSmRWko/VTwSXErJ0KmtLnVZmOBUrW9uWlMq3rlIz3Fa5rtfWq2owVPfncuWNDYCKSVySS4nAlukisi51XS/1eqFFWrV9GbTwr1LJrW6tqz6ocqttbuB7Apc6sUuNtpbbfURiHUoEVp3MpU3qUvqrRmK1JcCUvqqbHZTIfbUWSWheVIgpiz5A98P7Kl6g/v2B994DavBHPxHRuS7giOzEiRO9JmPpy9atnltmHj16FP369cNtt92G++67r8r2Bw8ejOuuuw4dO3bEwIED8fXXX2Pv3r346quvfNafMGECCgsLtZKdnR3oUyKqn/r0AbZsATp1AgA0KC7DoPE7MODJ/8F8yhXiztGAX7bji+kvVQxiJUnere3LLzmIJQoVRmTDTsAR2QcffBBDhgypsk7Lli21n48ePYo+ffqgR48emD17dsAdTEhIQEpKCvbt2+fzfrvd7jNSS0QAWreWv6IeNUqO8gHovPQIUrb8hW+e74DjPeWNDyKVPVGtyor8DZRjddmuRlb5Vt1A4aRN/j/ntFlQZpU/RtQcUXXpKq/IrE5FJFbNjVVzZj0jseoxKuXGak0qkVihhXn1y2up1dTltNTcV92tkhNbHgHlWLl1VERkXXbl+dmVHFglEiupkViT7jeV+sTLlEhssfyCtMguwpNff4KM3b9U1E1MBN5/H7jyShARUc0FPJCNjY1FbGxsjeoeOXIEffr0QVpaGubPnw/TGWynWFBQgOzsbCQkJAR8LRFBnjS0cCFw7bXA6NHAqVNodPg0hgzfhi23p+DbMW3RvMkJuaoygHUrX9ZYlZlZEcoCpw2UdWYjbHJKwkmrC+Xq5CaLMqg06wayJs9bv4TnZK6KdWN1t4CPbcPUx1DSHLRduJRjfwNXh24AqwxaKx+7HfLzcysDWdjUFALfA1jhUp5oifwCmE/Kx44Cgdu3b8I/v/0SDcsqrec7eLC80kTTpj6fExEFjwSvP4drrV0yhmGTvY4ePYorr7wSSUlJmDJlCv7880/k5uYiNzfXo17btm2xdOlSAMDJkycxfvx4bNy4EQcPHsSaNWswcOBAxMbG4iZOeiA6c5IkrzP788/AZZfJ54RAt8UHMfrGNWj15Z+eA0WqVR0OH8Z7C9/E0998XDGIbdYMWLRILhzEEhGdEcMme61YsQL79+/H/v370aJFC4/7RKVfmHv27EFhobwkkNlsxi+//IKFCxfi+PHjSEhIQJ8+fbB48WJERUUZ1VWic0fr1sCaNfIyXU88ARQXIzqvBFf+cy+SFh7DukcvwNHOjWFXluNyKrOhypTwqkn5Cl9dpstidmuRSXW5KS0Sq7/VRWbVyK3257S2y5bkcSz5ir4qkVctCmxVJnepEVi7fOxyKMuHqbcRcltlDdSIrNycS86wqDiOEMr1AkJJJZD8RGKFMplLKMtqQUkhsBaa0fTECYz/4mvc/MuPMFVOkhsxAnjlFebCEtU1RuWzMk5gGMMGsiNGjMCIESOqrVd5UBsREYHly5cb1SUiAuRtbceMAW68EXjgAeDrrwEArX/+E63v/BO/XR2HHx9ohT/b8o/HM9XwdAnuX7ke93y/BpGVdzO88ELgrbeYC0tUR3FDhPATtOW3iKiOSU0Fli2TB7JjxwJ79gAA2n77B9p++we2XJWCVcPb4VCHWG0ZrhJlnSoh1GipqIhMKpFXt8Vz0pdbO++55JUoVe9X2lLuNymbLqi5ssJzjwZAkrRcWC0Ca1MjsLpIrBqBVW7LdZHY8kgo9ZXbBko+rJIXC7sbZqsy6UxbVkt5fkourBaJPWmBo7QUd6/ciH+sXYXGxacq+hwVJW9S8dBDAHcqJCKqNRzIEp3r+vcHrrlG3kThmWcAJY+923eH0O27Q9h7SXOsvK09tl2ZHOKO1l3Rp4px9zeb8bc169H0ZKUBrMUirxjx5JNAXFzoOkhENcPUgrDDgSwRAVarPOAaPhyYPRuYPBn44w8AwIU/5eHCn/LwR/MofNWvE77u2wGnGstRxbJysxad1b4703JglVuL5626+YDJqm4UYPK4XMuVVSK0ul1zIUwVubHqqgRaBNbuGYnVbtUcWGU1Au24gRJdjVAeJEKOvpps8q3F6oJJtyqBS+mvu8yEFn/+hbu/3oyhP2xCVInTs6N33AFMmiTnJRMRkSE4kCWiChERwCOPAPffD7z7LvDqq8BvvwEA4vJO4J6FG3DPwg34ud15+Pqqjvis20XItZxbX5Wb3G5cvmMfhn69BVf+by9MlVd7MJmAIUOAxx7TNqIgojDD6GlY4UCWiLw5HMDIkcB99wGrVgHTpsn5tG45ctl59xF03n0E/zSvxKa2rfDtxe2w9sJ2yGweW2lVAvlW2wZW2WRA3apVi7RKnuvPmmxKA+qGCeoysurmBhapYmMDe9XrwuojsC4lAqvmwqKBHHk12+WNIOx2ebUGm5IXaza55YizEEja/xeuW78T163ZifMKjnu+XjabvBLBo48yAktEFEQcyBKRf5Ikb6Rw7bXA4cPyxgoffgjs3AkAsLrc6P3rfvT+dT+AL5DdtDHWtWmL79tciG3ntcaxyIZwKwFbl7orrjowNXlO/jKVq6kGykO7dRsNaLtzAW6rOoCFx62/gatbGbBKDeQBqyNCHrA2cCgbPdiUjR4s8nmTJAAhkLz/L1y6OhPXrtuNlof/8n59kpPllIx77mEOLFE9wFULwg8HskRUMy1aAI8/DkyYIG+s8OGHwEcfAQcPalWSCo7hzg0bceeGjQCAfc3j8EtiMn5JTMa+xon4rXkiSlFHt5QWAglHjqPTL0fQ+efD6LI1G83yT3rXM5uBvn3lAeyAAfIxERGFBAeyRBQYSQIuvlguL7wgR2eXLQNWrgTWrQPKyrSqF+T9gQvy/sDNO7YAAFyShINNmmFPXCJ+bxKPA02aIzu6GbIbNcVpi7wGVkXkQt0gQT6qnKrgUqK8bi0Sq5u0pURebZFyX6IayLtpNXKclm9tp9E4vxhtDuQiee9fOP+3PLTe8Sca/Xna//O+/HLg1luB229n9JWovuKqBWGHA1kiOnOSJE9q6tQJ+Ne/gJMngbVrge++A9avB7ZvB8rLtepmIdC6IA+tC/K8miqyRyA3qhHyImNQ0CAKfzVoiOOOSByPaIDjjkgUNYjAKZsDJxvYUBxhRanFgjKXCS6TCW7hhtnthsVZBntZGRqiGFGnS9Cs7AQaF51CYnEhmh47ifOOFyIutwiJOcfRoLjMqw8eHA7giiuAG24ABg0CEhNr+cUjorqGqQXhhwNZIqo9DRsC110nFwA4fVpOQ/jpJ3lQ+9NPcgS38m5XimjnaUQ7T+PC/Jwgd1rRsCGQni4PXi+7TP7Z4QhNX4iIqEY4kCUi40REyAPC9PSKc2VlwL59wK+/yrf79gGZmcChQ8DRoz4HubXKYgFSUoALLgDat5ejyZdcIv9s4Uci0TmNqQVhh5/aRBRcVqs8aGzf3vs+IYD8fCAvTy4FBfLxsWNyOXFCTl8oKZGLyyUvCSZJ8qQru12OokZFyaVpUyA2FmjeXM5rPe88ID6eA1YionqCn+ZEVHdIEtCsmVw6dAh1b4joHMMc2fBjCnUHiIiIiIjOBCOyRERERABzZMMQI7JEREREFJYYkSUiIiICGJENQxzIEhEREYGTvcIRUwuIiIiIKCwxIktEREQEMLUgDDEiS0RERERhiRFZIiIiIgCSEJBE7YdPjWiTZIzIEhEREdUh69atw8CBA5GYmAhJkvDpp59We83atWuRlpYGh8OBVq1aYdasWV51lixZgvbt28Nut6N9+/ZYunSpV50ZM2YgNTUVDocDaWlpWL9+fW08JcNwIEtEREQEVOTIGlECcOrUKXTu3BlvvvlmjepnZmZiwIAB6N27N7Zv347HH38cDz/8MJYsWaLV2bhxIwYPHoxhw4bh559/xrBhw3D77bfjxx9/1OosXrwYY8aMwRNPPIHt27ejd+/e6N+/P7KysgJ7AkEkCVG/4t1FRUWIiYlBYWEhoqOjQ90dIiIiqoFQ/v5WH/viu/4Ds81R6+27Skuw470nzui5SZKEpUuXYtCgQX7r/Otf/8Lnn3+O3bt3a+dGjRqFn3/+GRs3bgQADB48GEVFRfj666+1Ov369UPjxo3x4YcfAgC6d++OSy65BDNnztTqtGvXDoMGDcLkyZMD6newMCJLREREhIp1ZI0ogDxgrlycTmet9Hvjxo3IyMjwONe3b19s3boVZWVlVdbZsGEDAKC0tBTbtm3zqpORkaHVqYs4kCUiIiIKgqSkJMTExGiltqKcubm5iIuL8zgXFxeH8vJy5OfnV1knNzcXAJCfnw+Xy1VlnbqIqxYQERERAYavI5udne2RWmC322vtISRJ8nxIJXO08nlfdfTnalKnLuFAloiIiAjGb1EbHR1tSP5vfHy8V9Q0Ly8PFosFTZs2rbKOGoGNjY2F2Wyusk5dxNQCIiIiojDWo0cPrFy50uPcihUr0LVrV1it1irr9OzZEwBgs9mQlpbmVWflypVanbqIEVkiIiIioM5sUXvy5Ens379fO87MzMSOHTvQpEkTJCcnY8KECThy5AgWLlwIQF6h4M0338S4ceMwcuRIbNy4EXPnztVWIwCARx55BJdffjlefPFF3Hjjjfjss8+watUqfP/991qdcePGYdiwYejatSt69OiB2bNnIysrC6NGjTq7528gDmSJiIiI6pCtW7eiT58+2vG4ceMAAMOHD8eCBQuQk5PjsbZramoqli1bhrFjx2L69OlITEzEtGnTcMstt2h1evbsiUWLFuHJJ5/EU089hdatW2Px4sXo3r27Vmfw4MEoKCjApEmTkJOTg44dO2LZsmVISUkJwrM+M1xHloiIiEKuLqwjmzbYuHVkty0+s3VkqWrMkSUiIiKisMTUAiIiIiKgzuTIUs0xIktEREREYYkRWSIiIiKFEevIknE4kCUiIiICACHkYkS7ZAimFhARERFRWGJEloiIiAjGb1FLtY8RWSIiIiIKS4zIEhEREQFcfisMMSJLRERERGGJEVkiIiIiAJJbLka0S8ZgRJaIiIiIwhIjskREREQAc2TDEAeyRERERODyW+HI0NSCli1bQpIkj/LYY49VeY0QAhMnTkRiYiIiIiJw5ZVX4tdffzWym0REREQUhgzPkZ00aRJycnK08uSTT1ZZ/6WXXsKrr76KN998E1u2bEF8fDyuvfZanDhxwuiuEhER0blM3aLWiEKGMHwgGxUVhfj4eK00bNjQb10hBKZOnYonnngCN998Mzp27Ih33nkHxcXF+OCDD4zuKhERERGFEcMHsi+++CKaNm2Kiy++GP/5z39QWlrqt25mZiZyc3ORkZGhnbPb7bjiiiuwYcMGn9c4nU4UFRV5FCIiIqJAqTmyRhQyhqGTvR555BFccsklaNy4MTZv3owJEyYgMzMTb7/9ts/6ubm5AIC4uDiP83FxcTh06JDPayZPnoxnnnmmdjtORERERHVewBHZiRMnek3g0petW7cCAMaOHYsrrrgCF110Ee677z7MmjULc+fORUFBQZWPIUmSx7EQwuucasKECSgsLNRKdnZ2oE+JiIiIqGL5LSMKGSLgiOyDDz6IIUOGVFmnZcuWPs+np6cDAPbv34+mTZt63R8fHw9AjswmJCRo5/Py8ryitCq73Q673V6TrhMRERFRPRLwQDY2NhaxsbFn9GDbt28HAI9BamWpqamIj4/HypUr0aVLFwBAaWkp1q5dixdffPGMHpOIiIioJriObPgxbLLXxo0b8dprr2HHjh3IzMzERx99hPvvvx833HADkpOTtXpt27bF0qVLAcgpBWPGjMHzzz+PpUuXYufOnRgxYgQaNGiAO+64w6iuEhEREXH5rTBk2GQvu92OxYsX45lnnoHT6URKSgpGjhyJRx991KPenj17UFhYqB0/+uijOH36NEaPHo1jx46he/fuWLFiBaKioozqKhERERGFIUmI+vVnQlFREWJiYlBYWIjo6OhQd4eIiIhqIJS/v9XH7tF/EixWR623X15Wgo1f/5tjEwMYvo4sEREREZERDF1HloiIiChsGLVUVr367rtuYUSWiIiIiMISI7JERERE4PJb4YgRWSIiIiIKS4zIEhEREQGAW8jFiHbJEBzIEhEREQGc7BWGmFpARERERGGJEVkiIiIiABIMmuxV+02SghFZIiIiIgpLjMgSERERAYAQcjGiXTIEI7JEREREFJYYkSUiIiICN0QIR4zIEhEREVFYYkSWiIiICOA6smGIA1kiIiIiAJIQkAyYmGVEmyRjagERERERhSVGZImIiIgAwK0UI9olQzAiS0RERERhiRFZIiIiIjBHNhwxIktEREREYYkRWSIiIiKAy2+FIUZkiYiIiCgsMSJLREREBABCyMWIdskQHMgSERERAZCEXIxol4zB1AIiIiIiCkscyBIREREBFakFRpQAzZgxA6mpqXA4HEhLS8P69eurrD99+nS0a9cOERERaNOmDRYuXOhx/5VXXglJkrzKddddp9WZOHGi1/3x8fEB9z2YmFpAREREVIcsXrwYY8aMwYwZM9CrVy+89dZb6N+/P3bt2oXk5GSv+jNnzsSECRMwZ84cdOvWDZs3b8bIkSPRuHFjDBw4EADwySefoLS0VLumoKAAnTt3xm233ebRVocOHbBq1Srt2Gw2G/QsawcHskREREQAJLdcjGg3EK+++iruvfde3HfffQCAqVOnYvny5Zg5cyYmT57sVf/dd9/F/fffj8GDBwMAWrVqhU2bNuHFF1/UBrJNmjTxuGbRokVo0KCB10DWYrHU+ShsZUwtICIiIgqCoqIij+J0Or3qlJaWYtu2bcjIyPA4n5GRgQ0bNvhs1+l0wuFweJyLiIjA5s2bUVZW5vOauXPnYsiQIYiMjPQ4v2/fPiQmJiI1NRVDhgzBgQMHAnmKQceBLBERERFgeI5sUlISYmJitOIrupqfnw+Xy4W4uDiP83FxccjNzfXZ7b59++Ltt9/Gtm3bIITA1q1bMW/ePJSVlSE/P9+r/ubNm7Fz504t4qvq3r07Fi5ciOXLl2POnDnIzc1Fz549UVBQcKavqOGYWkBEREQUBNnZ2YiOjtaO7Xa737qSJHkcCyG8zqmeeuop5ObmIj09HUIIxMXFYcSIEXjppZd85rjOnTsXHTt2xKWXXupxvn///trPnTp1Qo8ePdC6dWu88847GDduXI2eY7AxIktEREQEVGxRa0QBEB0d7VF8DWRjY2NhNpu9oq95eXleUVpVREQE5s2bh+LiYhw8eBBZWVlo2bIloqKiEBsb61G3uLgYixYt8orG+hIZGYlOnTph37591dYNFQ5kiYiIiABIQhhWaspmsyEtLQ0rV670OL9y5Ur07NmzymutVitatGgBs9mMRYsW4frrr4fJ5DnU++ijj+B0OnHXXXdV2xen04ndu3cjISGhxv0PNqYWEBEREdUh48aNw7Bhw9C1a1f06NEDs2fPRlZWFkaNGgUAmDBhAo4cOaKtFbt3715s3rwZ3bt3x7Fjx/Dqq69i586deOedd7zanjt3LgYNGoSmTZt63Td+/HgMHDgQycnJyMvLw3PPPYeioiIMHz7c2Cd8FjiQJSIiIgLOePOCGrUbgMGDB6OgoACTJk1CTk4OOnbsiGXLliElJQUAkJOTg6ysLK2+y+XCK6+8gj179sBqtaJPnz7YsGEDWrZs6dHu3r178f3332PFihU+H/fw4cMYOnQo8vPz0axZM6Snp2PTpk3a49ZFkhBGvGOhU1RUhJiYGBQWFnokVBMREVHdFcrf3+pj90mbAIvFUf0FASovL8HqbZM5NjEAI7JEREREgDwpy4ANEVCvQoZ1Cyd7EREREVFYYkSWiIiICAh4hYFA2iVjMCJLRERERGGJEVkiIiIiQNm8wIhVC2q/SZIxIktEREREYYkRWSIiIiKgzqwjSzXHgSwRERERIC+9JRnULhmCqQVEREREFJYYkSUiIiICl98KR4zIEhEREVFYMmwgu2bNGkiS5LNs2bLF73UjRozwqp+enm5UN4mIiIhk6mQvIwoZwrDUgp49eyInJ8fj3FNPPYVVq1aha9euVV7br18/zJ8/Xzu22WyG9JGIiIiIwpdhA1mbzYb4+HjtuKysDJ9//jkefPBBSFLVUwLtdrvHtURERESG4/JbYSdoObKff/458vPzMWLEiGrrrlmzBs2bN8eFF16IkSNHIi8vz29dp9OJoqIij0JERERE9V/QBrJz585F3759kZSUVGW9/v374/3338d3332HV155BVu2bMFVV10Fp9Pps/7kyZMRExOjleraJyIiIvKJObJhJ+CB7MSJE/1O4lLL1q1bPa45fPgwli9fjnvvvbfa9gcPHozrrrsOHTt2xMCBA/H1119j7969+Oqrr3zWnzBhAgoLC7WSnZ0d6FMiIiIikjcuMKqQIQLOkX3wwQcxZMiQKuu0bNnS43j+/Plo2rQpbrjhhkAfDgkJCUhJScG+fft83m+322G32wNul4iIiIjCW8AD2djYWMTGxta4vhAC8+fPx9133w2r1Rrow6GgoADZ2dlISEgI+FoiIiKimuKGCOHH8BzZ7777DpmZmX7TCtq2bYulS5cCAE6ePInx48dj48aNOHjwINasWYOBAwciNjYWN910k9FdJSIiIqIwYvgWtXPnzkXPnj3Rrl07n/fv2bMHhYWFAACz2YxffvkFCxcuxPHjx5GQkIA+ffpg8eLFiIqKMrqrREREdC7j8lthx/CB7AcffFDl/aLSmxsREYHly5cb3SUiIiIiqgcMH8gSERERhQW3ACQDoqduRmSNErR1ZImIiIiIahMjskREREQAc2TDEAeyRERERAAAo3bh4kDWKEwtICIiIqKwxIgsEREREcDUgjDEiCwRERERhSVGZImIiIgAZZksLr8VThiRJSIiIqKwxIgsEREREQAIt1yMaJcMwYgsEREREYUlRmSJiIiIAK5aEIY4kCUiIiICONkrDDG1gIiIiIjCEiOyRERERABTC8IQI7JEREREFJYYkSUiIiIC5PRYQyKytd8kyRiRJSIiIqKwxIgsEREREcAc2TDEiCwRERERhSVGZImIiIgAwO0GYMB2sm5uUWsUDmSJiIiIAKYWhCGmFhARERFRWGJEloiIiAhgRDYMMSJLRERERGGJEVkiIiIiAHALGLJ7gZsRWaMwIktEREREYYkRWSIiIiIAQrghRO0vlWVEmyRjRJaIiIiojpkxYwZSU1PhcDiQlpaG9evXV1l/+vTpaNeuHSIiItCmTRssXLjQ4/4FCxZAkiSvUlJSclaPG2qMyBIREREB8uoCRuSzBrhqweLFizFmzBjMmDEDvXr1wltvvYX+/ftj165dSE5O9qo/c+ZMTJgwAXPmzEG3bt2wefNmjBw5Eo0bN8bAgQO1etHR0dizZ4/HtQ6H44wfty6QhKhfa0IUFRUhJiYGhYWFiI6ODnV3iIiIqAZC+ftbfeyrY4bBItlqvf1yUYpvC9+t8XPr3r07LrnkEsycOVM7165dOwwaNAiTJ0/2qt+zZ0/06tULL7/8snZuzJgx2Lp1K77//nsAckR2zJgxOH78eK09bl3A1AIiIiKiICgqKvIoTqfTq05paSm2bduGjIwMj/MZGRnYsGGDz3adTqdHZBUAIiIisHnzZpSVlWnnTp48iZSUFLRo0QLXX389tm/fflaPWxdwIEtEREQEAG63cQVAUlISYmJitOIrypmfnw+Xy4W4uDiP83FxccjNzfXZ7b59++Ltt9/Gtm3bIITA1q1bMW/ePJSVlSE/Px8A0LZtWyxYsACff/45PvzwQzgcDvTq1Qv79u0748etC5gjS0RERBQE2dnZHqkFdrvdb11JkjyOhRBe51RPPfUUcnNzkZ6eDiEE4uLiMGLECLz00kswm80AgPT0dKSnp2vX9OrVC5dccgneeOMNTJs27Ywety5gRJaIiIgIqNii1ogCebJV5eJrIBsbGwuz2ewVBc3Ly/OKlqoiIiIwb948FBcX4+DBg8jKykLLli0RFRWF2NhYn9eYTCZ069ZNi8ieyePWBRzIEhEREdURNpsNaWlpWLlypcf5lStXomfPnlVea7Va0aJFC5jNZixatAjXX389TCbfQz0hBHbs2IGEhISzftxQYmoBEREREQDhdkNIod8QYdy4cRg2bBi6du2KHj16YPbs2cjKysKoUaMAABMmTMCRI0e0tWL37t2LzZs3o3v37jh27BheffVV7Ny5E++8847W5jPPPIP09HRccMEFKCoqwrRp07Bjxw5Mnz69xo9bF3EgS0RERFSHDB48GAUFBZg0aRJycnLQsWNHLFu2DCkpKQCAnJwcZGVlafVdLhdeeeUV7NmzB1arFX369MGGDRvQsmVLrc7x48fx97//Hbm5uYiJiUGXLl2wbt06XHrppTV+3LqI68gSERFRyNWFdWSvihhs2Dqy351ezLGJARiRJSIiIgLkXb2k0O/sRTXHyV5EREREFJYYkSUiIiIClMhp7U/2YkTWOIzIEhEREVFYYkSWiIiICIBwCwgDcmTr2bz6OoURWSIiIiIKS4zIEhEREQGAcMOYHFkD2iQAjMgSERERUZhiRJaIiIgIzJENRxzIEhEREQFMLQhD9W4gq/7VU1RUFOKeEBERUU2pv7dDGb0sRxlgwMOXo6z2GyUA9XAge+LECQBAUlJSiHtCREREgTpx4gRiYmKC+pg2mw3x8fH4PneZYY8RHx8Pm81mWPvnKknUs8QNt9uNo0ePIioqCpIkhbo7higqKkJSUhKys7MRHR0d6u4E3bn+/AG+Buf68wf4Gpzrzx+of6+BEAInTpxAYmIiTKbgz0UvKSlBaWmpYe3bbDY4HA7D2j9X1buIrMlkQosWLULdjaCIjo6uFx9eZ+pcf/4AX4Nz/fkDfA3O9ecP1K/XINiR2MocDgcHmmGIy28RERERUVjiQJaIiIiIwhIHsmHIbrfj6aefht1uD3VXQuJcf/4AX4Nz/fkDfA3O9ecP8DUgAurhZC8iIiIiOjcwIktEREREYYkDWSIiIiIKSxzIEhEREVFY4kCWiIiIiMISB7JEREREFJY4kK3j1qxZA0mSfJYtW7b4vW7EiBFe9dPT04PY89rVsmVLr+fz2GOPVXmNEAITJ05EYmIiIiIicOWVV+LXX38NUo9rz8GDB3HvvfciNTUVERERaN26NZ5++ulqt1IM938DM2bMQGpqKhwOB9LS0rB+/foq669duxZpaWlwOBxo1aoVZs2aFaSe1r7JkyejW7duiIqKQvPmzTFo0CDs2bOnymv8fVb89ttvQep17Zk4caLX84iPj6/ymvr0/gO+P/MkScIDDzzgs359ev+JAlHvtqitb3r27ImcnByPc0899RRWrVqFrl27Vnltv379MH/+fO3YZrMZ0sdgmTRpEkaOHKkdN2zYsMr6L730El599VUsWLAAF154IZ577jlce+212LNnD6Kioozubq357bff4Ha78dZbb+H888/Hzp07MXLkSJw6dQpTpkyp8tpw/TewePFijBkzBjNmzECvXr3w1ltvoX///ti1axeSk5O96mdmZmLAgAEYOXIk3nvvPfzwww8YPXo0mjVrhltuuSUEz+DsrF27Fg888AC6deuG8vJyPPHEE8jIyMCuXbsQGRlZ5bV79uzx2K60WbNmRnfXEB06dMCqVau0Y7PZ7LdufXv/AWDLli1wuVza8c6dO3Httdfitttuq/K6+vL+E9WYoLBSWloqmjdvLiZNmlRlveHDh4sbb7wxOJ0KgpSUFPHaa6/VuL7b7Rbx8fHihRde0M6VlJSImJgYMWvWLAN6GFwvvfSSSE1NrbJOOP8buPTSS8WoUaM8zrVt21Y89thjPus/+uijom3bth7n7r//fpGenm5YH4MpLy9PABBr1671W2f16tUCgDh27FjwOmaQp59+WnTu3LnG9ev7+y+EEI888oho3bq1cLvdPu+vT+8/USCYWhBmPv/8c+Tn52PEiBHV1l2zZg2aN2+OCy+8ECNHjkReXp7xHTTQiy++iKZNm+Liiy/Gf/7znyq/Ws/MzERubi4yMjK0c3a7HVdccQU2bNgQjO4aqrCwEE2aNKm2Xjj+GygtLcW2bds83jsAyMjI8Pvebdy40at+3759sXXrVpSVlRnW12ApLCwEgBq95126dEFCQgKuvvpqrF692uiuGWbfvn1ITExEamoqhgwZggMHDvitW9/f/9LSUrz33nu45557IElSlXXry/tPVFMcyIaZuXPnom/fvkhKSqqyXv/+/fH+++/ju+++wyuvvIItW7bgqquugtPpDFJPa9cjjzyCRYsWYfXq1XjwwQcxdepUjB492m/93NxcAEBcXJzH+bi4OO2+cPX777/jjTfewKhRo6qsF67/BvLz8+FyuQJ673Jzc33WLy8vR35+vmF9DQYhBMaNG4fLLrsMHTt29FsvISEBs2fPxpIlS/DJJ5+gTZs2uPrqq7Fu3bog9rZ2dO/eHQsXLsTy5csxZ84c5ObmomfPnigoKPBZvz6//wDw6aef4vjx41UGMOrT+08UkFCHhM9VTz/9tABQZdmyZYvHNdnZ2cJkMomPP/444Mc7evSosFqtYsmSJbX1FM7ambwGqo8//lgAEPn5+T7v/+GHHwQAcfToUY/z9913n+jbt2+tP5czcSbP/8iRI+L8888X9957b8CPVxf/Dfhy5MgRAUBs2LDB4/xzzz0n2rRp4/OaCy64QDz//PMe577//nsBQOTk5BjW12AYPXq0SElJEdnZ2QFfe/3114uBAwca0KvgOnnypIiLixOvvPKKz/vr8/svhBAZGRni+uuvD/i6+vL+E1WFk71C5MEHH8SQIUOqrNOyZUuP4/nz56Np06a44YYbAn68hIQEpKSkYN++fQFfa5QzeQ1U6uz7/fv3o2nTpl73qzOcc3NzkZCQoJ3Py8vzityESqDP/+jRo+jTpw969OiB2bNnB/x4dfHfgC+xsbEwm81e0deq3rv4+Hif9S0Wi89/H+HioYcewueff45169ahRYsWAV+fnp6O9957z4CeBVdkZCQ6derk999ufX3/AeDQoUNYtWoVPvnkk4CvrS/vP1FVOJANkdjYWMTGxta4vhAC8+fPx9133w2r1Rrw4xUUFCA7O9tjUBdqgb4GlW3fvh0A/D6f1NRUxMfHY+XKlejSpQsAOc9s7dq1ePHFF8+sw7UskOd/5MgR9OnTB2lpaZg/fz5MpsCzgurivwFfbDYb0tLSsHLlStx0003a+ZUrV+LGG2/0eU2PHj3wxRdfeJxbsWIFunbtekb/X0JNCIGHHnoIS5cuxZo1a5CamnpG7Wzfvr3Ov9814XQ6sXv3bvTu3dvn/fXt/a9s/vz5aN68Oa677rqAr60v7z9RlUIdEqaaWbVqlQAgdu3a5fP+Nm3aiE8++UQIIcSJEyfEP//5T7FhwwaRmZkpVq9eLXr06CHOO+88UVRUFMxu14oNGzaIV199VWzfvl0cOHBALF68WCQmJoobbrjBo17l10AIIV544QURExMjPvnkE/HLL7+IoUOHioSEhLB7DdR0gquuukocPnxY5OTkaKWy+vRvYNGiRcJqtYq5c+eKXbt2iTFjxojIyEhx8OBBIYQQjz32mBg2bJhW/8CBA6JBgwZi7NixYteuXWLu3LnCarWeURpOXfCPf/xDxMTEiDVr1ni838XFxVod/Wvw2muviaVLl4q9e/eKnTt3iscee0wAqPOpJL7885//FGvWrBEHDhwQmzZtEtdff72Iioo6Z95/lcvlEsnJyeJf//qX1331+f0nCgQHsmFi6NChomfPnn7vByDmz58vhBCiuLhYZGRkiGbNmgmr1SqSk5PF8OHDRVZWVpB6W7u2bdsmunfvLmJiYoTD4RBt2rQRTz/9tDh16pRHvcqvgRDyElxPP/20iI+PF3a7XVx++eXil19+CXLvz978+fP95tBWVt/+DUyfPl2kpKQIm80mLrnkEo+lp4YPHy6uuOIKj/pr1qwRXbp0ETabTbRs2VLMnDkzyD2uPf7e78r/vvWvwYsvvihat24tHA6HaNy4sbjsssvEV199FfzO14LBgweLhIQEYbVaRWJiorj55pvFr7/+qt1f399/1fLlywUAsWfPHq/76vP7TxQISQghghwEJiIiIiI6a1x+i4iIiIjCEgeyRERERBSWOJAlIiIiorDEgSwRERERhSUOZImIiIgoLHEgS0RERERhiQNZIiIiIgpLHMgSERERUVjiQJaIiIiIwhIHskREREQUljiQJSIiIqKw9P/NxWW9n5/jmQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\"\"\" Plot the copula \"\"\"\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "\n", "# get the axis values for the copula\n", "xvals, yvals = kde.axes\n", "\n", "# find the 99.5th percentile probability contour of the original data\n", "pval = 1 - 0.995\n", "p_contour_val = fastkde.plot.calculate_probability_contour(\n", " pdf = kde.pdf[..., np.newaxis],\n", " axes = [[0]] + kde.axes,\n", " pvals = pval)\n", "\n", "# mask the copula wherever the PDF is less than the 99.95th percentile\n", "mask = kde.pdf < p_contour_val\n", "copula = np.ma.masked_where(mask, copula)\n", "\n", "# plot the copula\n", "cplt = ax.pcolormesh(xvals, yvals, copula, cmap='viridis', shading='auto')\n", "\n", "# add a colorbar\n", "cbar = plt.colorbar(cplt, ax=ax)\n", "cbar.set_label('Copula Density', rotation=270, labelpad=20)\n", "\n", "# plot the 99.5th percentile contour\n", "ax.contour(kde.axes[0], kde.axes[1], kde.pdf, levels=p_contour_val, colors='red', linewidths=2)\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "8c5d5955", "metadata": {}, "source": [ "```{note}\n", "The copula density is estimated by first estimating the marginal densities of each variable, and then estimating the joint density of the data in the transformed space. The copula density is then obtained by dividing the joint density by the product of the marginal densities.\n", "\n", "This leads to large variations in the copula density in regions where the marginal densities are small; this is especially true in regions where data are sparse. Note in the above that the copula density is masked such that it is only shown in the region where 99.95% of the data lies (as estimated by the `calculate_probability_contour` function). \n", "```\n", "\n", "Within the region where the copula density is shown, it is approximately constant, with values near 1.0. But there are large variations along the edge. This illustrates the difficulty in empirically estimating the copula density.\n" ] } ], "metadata": { "kernelspec": { "display_name": "fastkde_dev", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.10" } }, "nbformat": 4, "nbformat_minor": 5 }